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ABSTRACT 


With the increasing opportunities for research in a 
microgravity environment, there arises a need for understanding 
fluid mechanics under such conditions. In particular, a number 
of material processing configurations involve fluid-fluid 
interfaces which may experience instabilities in the presence 
of external forcing. In a microgravity environment, these 
accelerations may be periodic or impulse-type in nature. This 
research investigates the behavior of a multi-layer idealized 
fluid configuration which is infinite in extent. The analysis 
is linear, and each fluid region is considered inviscid, 
incompressible, and immiscible. 

An initial parametric study of configuration stability in 
the presence of a constant acceleration field is performed 
The zero mean gravity limit case serves as the base state for 
the subsequent time-dependent forcing cases. A stability 
analysis of the multi-layer fluid system in the presence of 
periodic forcing is investigated. Floquet theory is utilized. 
K parameter study is performed, and regions of stability are 
identified. For the impulse-type forcing case, asymptotic 
stability is established for the configuration. Using numerical 
integration, the time response of the interfaces is determined. 
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CHAPTER 1 


INTRODUCTION 


1.1 Literature review 

The microgravity environment aboard the space shuttle has 
given rise to an number of research opportunities which will 
increase when space station becomes operational. In 
particular, materials processing , which generally involves 
fluid configurations, will involve processes which exhibit 
significantly different behavior in a microgravity environment. 
The gravity- induced fluid— thermal flows; ie. buoyancy-driven 
convection in liquids, will no longer contribute. This 
physical phenomenon masked the presence of thermo-capillary 
flows, which will now assume a greater role in a microgravity 
environment 1 . 

The effect of gravity has been greatly reduced in 
low-gravity aircraft flights and drop tubes which provide short 
periods of microgravity, sufficient for some research, but 
certainly too brief for most materials processing experiments. 
The advent of extended spaceflight has dramatically increased 
the opportunities for long-duration research and development in 
space. There are numerous technological applications which are 
envisioned in a microgravity environment. 

The growth of crystals for electronic materials has not 
reached theoretical performance limits due to defects caused in 
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part by the presence of gravity. During the spacelab missions, 
scientists were able to monitor growth of a crystal through 
each stage of its formation. In earth-grown crystals, it can 
be observed where the seed crystal stops and the new growth 
begins. The introduction of such a defect was not detected in 
space due to the lack of gravity-induced convection 20 . 

The great reduction in convection is also relevant to 
metallurgical manufacturing. A microgravity environment 
provides greater understanding of how liquefied metals diffuse 
through each other prior to solidification. Such knowledge is 
important for the production of improved and novel alloys. 

Containerless processing makes possible the production of 
much improved glasses and ceramics. In such a process, the 
sample is suspended and manipulated by acoustical and 
electromagnetic forces without the contamination of a 
container. Large samples can be dealt with in a microgravity 
environment 20 . 

Biological processing also benefits from space. Large, 
pure crystals allow analysis of many unknown protein structures 
which are essential to the design of new and improved drugs. 
There is also effort towards the separation and purification of 
biological substances for pharmaceutical purposes 17,21 . 

In the absence of gravity, fluid behavior which might 
normally be hidden by gravity-driven flows in a terrestrial 
environment can be observed and analyzed. Drop and fluid 
column dynamics in microgravity permit experimentation of basic 
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fluid physics theories. Experiments have been performed 
concerning the stability of liquid bridges in a short term 
microgravity on a rocket 18 . Also, experimental work has been 
done on board Spacelab to investigate the shapes of rotating 
free drops in a microgravity environment 28 . in fact, the fluid 
configurations of drop dynamics and liquid columns will occur 
not only in fundamental studies, but also in materials 
processing applications. For example, the proposed 
solidification of novel alloys could take place in an acoustic 
levitation chamber, with the liquid sample having a drop 
configuration. 

A float zone configuration can be utilized in the growth 
of crystals. The float zone itself is modeled by a liquid 
column. In materials processing applications, heat and mass 
transfer effects are present in addition to the fluid dynamics. 
In fact, Marangoni convection would occur in the liquid column 
in a realistic processing scenario. It is currently thought 
that this convection could be reduced via the addition of a 
surrounding layer of fluid around the float zone. This would 
result in a multi-layer, compound fluid column configuration 2 . 

The environment on board a spacecraft is not strictly a 
microgravity environment. Rather, residual accelerations exist 
which could affect any ongoing materials science or space 
processing experiments. A recent summary 22 indicates that 
on board the space shuttle, accelerations include those in the 
frequency range up to ten hertz, with acceleration levels 
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In addition to 


ranging from 10 ^*gearth to 10 *g«»rth. 
periodic accelerations (g-jitter), residual accelerations may 
be of impulse type, due to such causes as station-keeping 
maneuvers and astronaut motion 27 . 

Most processes involve fluid dynamics, and in particular, 
fluid interfaces. This study does not investigate a specific 
process per se, but instead considers the stability of 
initially planar fluid interfaces. 

Previous work on fluid interfaces in microgravity has 
focused predominantly on the application of fuel slosh in 
tanks. Most recently, this has included work done by Hung et 
al 9 , which considered g-jitter in a slosh tank. A brief 

review of earlier work, as well as an extension of the previous 
efforts, was given by Gu et al 7 . These investigations all 
involved liquid in a container of specified shape with a free 
surface. 

The stability of a single planar free surface subjected to 
periodic forcing in the direction perpendicular to the 
interface has been investigated 3,6 . Both studies were done in 
a 1-g ambient environment and required the use of a container. 
In the work of Benjamin and Ursell 3 , the container was 
cylindrical in shape. The analysis led to a Mathieu equation 
which governed the time— dependent amplitude of the disturbance. 
They were able to make statements concerning the interface 
stability based upon known mathematical properties of Mathieu 
equations. The case of a rectangular container has been 
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addressed recently by Gu 6 , and the results extended into the 
nonlinear regime. Both of these investigations utilized an 
inviscid analysis. 

Viscous effects on the stability properties have been 
investigated recently in idealized infinite or semi-infinite 
configurations which have one fluid-fluid interface 8,12 . The 
forcing was periodic and directed perpendicular to the 
interface. The work of Jacqmin and Duval 12 assumed a zero mean 
g-level and pertained to a microgravity environment. A Floquet 
analysis was applied to the fluid system for the case of 
sinusoidal forcing. Stability boundaries were obtained from 
the results. 

Recently, Jacqmin has studied the stability of an 
oscillated fluid with a uniform density gradient 11 . The case 
of forcing perpendicular to the density gradient was 
investigated. Such a problem involves the Kelvin-Helmholtz 
instability. 


1.2 Objectives 

This research will consider a multi-layer fluid 
configuration unbounded in space. Multi-layer fluid 
configurations are finding applications in materials processing 
scenarios in microgravity. Although the infinite, multiple 
layer fluid system is not physically realistic, it is the 
logical extension of work done previously in the one interface, 
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infinite system. As such, it will provide insight into the 
behavior of configurations with multiple interfaces in the 
presence of forcing. 

A layer of finite height will be situated between two 
semi-infinite layers of fluid. The analysis is linear, and 
each fluid region is considered inviscid, irrotat ional , 
incompressible, and immiscible. A normal mode approach in the 
spatial variables will be assumed. Surface tension is the only 
property of each interface and is taken to be constant. 
Jacqmin and Duval showed that the presence of even weak surface 
tension can overwhelm the effects of viscosity, making the 
viscous analysis of secondary importance 12 . 

The objective is to investigate the configuration behavior 
in the presence of microgravity environment acceleration 
fields. As stated previously, these may manifest as periodic 
or non-periodic impulse-type accelerations. It is recognized 
that in practice, the true acceleration field will be random in 
magnitude and orientation. Two subcases will be investigated: 
1) periodic forcing directed normal to the interface (a cosine 
forcing function will be assumed), 2) non-periodic but 
time-dependent normal impulse forcing. 

As a preliminary step to this investigation, the stability 
of the configuration in the presence of a constant acceleration 
field will be investigated (Chapter 2). Regions of stability 
and various parametric trends will be established. The zero 
mean gravity limit case will ultimately serve as the base state 
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for the investigation of the time-dependent acceleration cases. 

The periodic forcing case results in a system of four 
ordinary differential equations (in time) with periodic 
time-dependence. Such a problem is well-posed for application 
Floquet theory 5, 25 in which the time-dependent coefficients 
are expressed in terms of a Floquet exponent. Previous work in 
fluid mechanics has utilized Floquet theory 5,12 . It is, 
however, more generally applied to dynamical systems 25 . One 
recent application involved the analysis of a 
spin-stabilized satellite in orbit 26 . In the periodic forcing 
case, the system of fluid equations can be converted to an 
infinite algebraic eigensystem. The nature of the real 
component of the eigenvalue will determine configuration 
stability. The effect of six non-dimensional parameters will 
be investigated. 

For the non-periodic case, asymptotic stability will be 
established according to mathematical theory. The system 
results in four linear differential equations which will be 
integrated numerically. The time response of each interface 
will be determined, and parametric trends will be discerned. 
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CHAPTER 2 


MULTI-LAYER FLUID CONFIGURATION STABILITY IN THE PRESENCE 
OF CONSTANT ACCELERATION FIELDS 


2.1 Problem description 

Prior to an analysis of stability of a multi-layered 
configuration in the presence of time-dependent forcing, cases 
will be considered in which the body force is due entirely to a 
constant gravitational acceleration. The limit cases of 
l*g««rth and 0*g*«rth are studied, as well as various 
intermediate values. The 0*g»«rth mean state will ultimately 
serve as a basis for investigating the effects of residual 
accelerations in Chapters 3 and 4. 

The configuration to be considered is comprised of three 
horizontal fluid layers. No rigid boundaries are present. The 
layers extend to infinity in the horizontal directions. 
The top and bottom layers are considered to be semi-infinite 
in nature, while the middle layer has a finite height. The 
geometry of the figure is given by Figure 2.1. 

The base state is one of zero mean motion in each of the 
three fluid regions. The fluids are immiscible and will be 
taken as inviscid, irrotational and incompressible. Surface 
tension is a property of the interfaces. A normal mode 
approach is assumed. That is, the small amplitude 

perturbations are wavelike in nature. 

The fluids considered in this study are water, air, and 
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Configuration Geometry 


region 2 p 2 , to 


y = h + erin 



Yn 


region 1 pi, <|>i 


y = erjm 



Ym 


region 3 P3. to 


p = density of subscripted region 
$ = potential function of subscripted region 
Y = surface tension of subscripted interface 
"H = perturbation of subscripted interface 


Figure 2.1 


silicone oil. Four different cases are examined in the 
following configurations: 

CASE 1: air/silicone oil/water (region 2/region 1/region 3) 
CASE 2: air/water/air " 

CASE 3: air/silicone oil/air " 

CASE 4: water/silicone oil/water " 

The parameters to be varied include height of the middle 
slab, wave number of the interfacial perturbation, and the 
value of the constant gravitational acceleration. By varying 
these guantities, the propagation speed of the perturbations 
can be calculated for different parameter conditions. A 
positive value of the imaginary component of the propagation 
speed will indicate an instability on the fluid system. 

Several of the cases to be investigated have a 
configuration such that the density of the upper fluid is 
greater than that of the lower fluid, giving rise to a motion 
driven by gravity. This type of instability is known as the 
Rayleigh-Taylor instability. It will be shown that the growth 
rate of these instabilities is determined by the nature of the 
solution to the dispersion relation. More specifically, if a 
certain configuration generates non-zero imaginary components 
of the propagation speed, then depending on the sign of the 
quantity, the Rayleigh-Taylor instability will occur. 
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2.2 Equation development 
2 . 2a Governing equations 

As stated previously, a normal mode perturbation has been 
utilized. The waveform of the disturbance is given by the 
following: 


n(x,t) - Ce lk(X “Ct) (2.1) 

where <; = amplitude (small) 

k ■ wave number (real number) 
c » propagation speed (complex number) 

7) ■ interface shape 

The governing equations of incompressible fluid mechanics 
are the continuity equation and the momentum equation. The 
analysis is inviscid, irrotational , and linear. Linearization 
is done about a base state of zero mean motion. The following 
equations govern fluid behavior. 

V*u' - 0 (2.2) 


au' 

P — - -7p' (2.3a) 

at 


o 


“ 7 Pmean + 


(2.3b) 
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Note the momentum equation has been 
(2.3a) and mean (2.3b) equations 
omitted for perturbation values. 

A potential function, <p , 
Substitution into equation (2.2) 

7 2 <p =* 0 


split into the perturbation 
Henceforth, primes will be 

with u - V<p, is defined, 
yields Laplace's equation. 

(2.4) 


Laplace's equation must be solved in all three regions. 
Separation of variables yields the following solutions for the 
potential functions. 


- [ Ae ky ♦ B.' ky ] . lk < x - ct > 
* 2 - Ce- ky . 11 ' ( *- ct > 


- D« kY « ik(1< ' ct) 


(2.7) 


where A, B, C, and D are constant coefficients. 

Furthermore, at each interface the perturbation is defined 
as 
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ik (x-ct ) 


( 2 . 8 ) 


t? 


II 


Ee 


( 2 . 9 ) 


is obtained by applying three 
interface. These three conditions 
boundary condition, (2) the 
matching of the normal component of the velocity, and (3) the 
normal force balance. 

The kinematic condition states that a particle of fluid 
which is at some point on the surface will always remain on 
that surface. This can be written as: 


'III 


= Fe ik(X_Ct) 


2 . 2b Boundary conditions 
The dispersion relation 
boundary conditions at each 
are: (1) the kinematic 


p (y- CT » - o 

Dt 


By converting into Euler ian form, and noting that x, y, 
and t are independent and that the waveform depends solely on x 
and t, the equation (after linearization) becomes: 
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( 2 . 11 ) 


— (x,t) = — (x,y +en,t) (y = 0 or h) 

at ay 


By applying a Taylor expansion and again neglecting 
quadratically small terms, the kinematic boundary condition at 
each interface simplifies to: 


^ (x,h, t) - ( x , t ) 

ay at 


( 2 . 12 ) 


£4 (X;0 ,t) - (x,t) 

ay at 


( 2 . 13 ) 


Imposition of the condition that the normal component of 
velocity be continuous across the interface yields: 



3* 9*2 

at y - h 

( 2 . 14 ) 

— 

ay ay 




9*! 9*3 

at y - 0 

( 2 . 15 ) 


9y 9y 




Finally, the third 

boundary condition is imposed. 

Taking 
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into account the surface tension, the normal force balance 
takes the following form: 

^lower"* ^upper = (2.16) 

where 7 = surface tension 

ft = the outward pointing normal to the interface 


By noting that the perturbation is a function of x and t 
only, the final boundary conditions can be derived at each 
interface. 




+ 9 ( ‘W>ni 


a 2 „ 


“7. 


II 


11 ax 2 


(2.17) 




+ g(P 1 -P 3 )7 )lII 


-7 


III 


ax 


hi 

2 


(2.18) 


2.2c Dispersion relation 

These six algebraic equations (2.12-2.15,2.17,2.18) form a 
homogeneous system in unknowns A,B,C,D,E, and F. In order to 
have a solution, the determinant corresponding to this system 
must equal zero. This results in the following dispersion 
relation. 
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2k_h 


(P 1 +P 3 ) (P 1 +P 2 ) e + (p 2 -p 1 ) (P 1 -P 3 ) 


{<V P 3>( 1(p 2 -p 1 )-» I i it) + (P! + p 2 >( l ( Pi"p3 ) - T iu k )} e 

+ (P 2 ~P 1 ) ( ?( p 3 _p i> + 7 IH k ) + (p l" p 3 } ( 1 (P 2 ” P 1 ) " T II k ^ 


2kh 


{( lC p i _p 3) _7 m k ) ( l (p 2" p l )_r II^} 


2 Jch 


+ ( Sc (p 2~ p l* ~ 7 II k ) ( ?< p 3 _p i) _ 7 in k ) 


( 2 . 19 ) 


Thus, this propagation speed, c, is given by the solution 
of a fourth order polynomial. It is also the eigenvalue. 
From the theory of roots of a polynomial, it is readily seen 
that there exists the possibility of complex roots to the 
dispersion relation. The propagation speed can be written as a 
complex number. 


c 


c * + 1C . 


( 2 . 20 ) 


Hence, the perturbation equations can actually be written as: 
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( 2 . 21 ) 


'II 


T) 


III 


^ikx (-ike t) (kc t) 
e e R e i ' 

_ ikx (-ike t) (kc t) 
Fe e v r ' e ' i ' 


( 2 . 22 ) 


The first two exponential terms of each equation are 

oscillatory in nature. The third exponential factor is a real 
number. An imaginary component equal to zero implies a neutral 
disturbance, while if the value is less than zero, the 
exponential term decays in time, and the system remains stable. 
However, if this imaginary component, c , is positive, the 
exponential term grows in time, resulting in an instability. 
This case is known as the Rayleigh-Taylor instability. An 
analytical limit case can be obtained from the full dispersion 
relation for the special case in which the ratios of the top 
and bottom densities to the middle density are negligibly 
small . 

[.*«*-!]=«♦ [(.*»♦!) ('Isg - IlSI* )]c 2 

♦ [(•*“-») ( I - ) (1 - )] (*.») 
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For such cases, the configuration will remain stable if the 


following inequality holds true. 



The scope of this study is to analyze the four 

previously stated cases under various parameter conditions. 
That is, by allowing the parameters to vary over a specified 
range, the roots of the dispersion relation can be calculated 
;and hence, interface stability can be determined. The 
parameters that are considered are the height of the middle 
layer, the wave number, and the value of the gravitational 
constant. For our ultimate purposes, we are most interested in 
the case in which the time-independent gravitational body force 
is zero. 


2 . 3 Results 

The dispersion relation was solved numerically using the 
DZPORC routine of the IMSL library. The DZPORC routine makes 
use of the Jenkins-Traub three-stage algorithm 13 , in which the 
roots are computed one at a time for real roots and two at a 
time for complex conjugate pairs. As the roots are found, the 
real root or quadratic factor is removed by polynomial 
deflation. 
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The fourth order polynomial (in c) has four roots. Because 
of the nature of the dispersion relation, the roots 
were generated in pairs. That is, for any given solution, 
there exist two pairs of roots, where each pair consists of 
the positive and negative values of a number. Physically, for 
real roots, this means the perturbation may propagate in either 
the positive or negative direction. For imaginary roots, it 
implies an instability will occur since these roots occur in 
complex conjugate pairs. If all the roots are real, the system 
will be stable. 

Figure 2.2 shows the four roots of the dispersion relation 
for Case 1 (air/silicone oil/water) . The roots are plotted 
over a range of gravity ratios (g/g«»rth) from 0 to 1.0. As is 
expected, since heavier fluids underlie lighter ones this case 
is stable for all parametric conditions. (Note that the 
non-zero roots are exclusively real.) A less dense fluid above 
a more dense fluid is stable to small perturbations in the 
presence of constant gravitational forcing. 

Figures 2.3, 2.4, and 2.5 show the dispersion solution for 
Cases 2,3, and 4, respectively. Each of these cases reveals 
the presence of a positive imaginary root, which in turn, 
implies an unstable configuration. This behavior is expected 
as each case involves a more dense fluid above a less dense 
fluid in its configuration. 

Since an instability depends solely on the presence of 
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positive imaginary roots, the subsequent figures will display 
these particular roots exclusively. 

The effect of wave number, (k) , on configuration stability 
is elucidated in Figures 2. 6-2.8. As k values increase, the 
configuration becomes more stable. Since k is inversely 
proportional to wavelength, the configuration is unstable to 
long wavelength perturbations. The restoring force required to 
maintain stability is greater in the long wavelength regime. 

Note that all cases are stable at 0*g««rth. The results 
of Case 1 do not appear since the configuration is stable for 
all parameter space. 

The effect of surface tension can be readily seen by 
comparing Figure 2.6 with 2.7, where the middle layers are 
water and silicone oil, respectively. Thus, while their 
densities are effectively the same, the water-air interfaces 
have surface tension values nearly three times that of the 
oil-air interfaces. With the increased restoring force, it is 
expected and confirmed that Figure 2.6 will be more stable than 
Figure 2.7. In the water case (Fig. 2.6), for a k value of 3, 
the system is stable up to g*0 . 65*g««rth. For the oil case 
(Fig. 2.7), for k*3, the configuration becomes unstable at 
g*0 .23 *ge«rth • 

From Figure 2.9, it is tempting to conclude that the 
middle slab height, (h) , has no effect on the stability of the 
configuration. This conclusion is valid for values of h which 
are large in comparison to the wavelength (recall X » 2n/k) . 
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When the nondimensional quantity, h/A, is greater than or on 
the order of one, middle layer height has little effect on the 
stability. In Figure 2.9, this corresponds roughly to values 
of h e 5cm (for k^lcm -1 ) . For h-lcm, the quantity, h/A, equals 
0.16 which is less than 0(1). From Figure 2.9 it is seen that 
this height is associated with the fastest growing 
instabilities. 

The effect of middle layer height is even more dramatic in 
Figures 2.10 and 2.11. The fastest growing instabilities for a 
given wavelength perturbation are associated with 
configurations with the smallest values of h/A. In Figure 
2.10, the smallest value of h/A equals 0.04 (corresponding to 
h«0.25cm, k-l.Ocm” 1 ). Note that this value relates to the 
fastest growing instability. 

The limit case (eq. 2.24) simulates a liquid layer 
situated between two layers of a gas, and its accuracy can be 
verified by comparing it to either Figure 2.6 or 2.7. 
According to (2.24), for Case 2 (air/water/air), and h-lcm, the 
instability should originate at g/g*4rth»0. 073 for k-1, and at 
g/g«*rth “0. 661 for K— 3. 

For Case 2 (air/silicone oil/air) , and h-lcm, the 
instability should start at g/g#4rth-0.027 for k-1, 
g/g««rth-0. 2 30 for k-3, and g/g«.rth-0 . 657 for k-5. The 
results from Figures 2.6 and 2.7 confirm these values. 

It is seen that in the case of zero gravity, each 
configuration remains stable. Although we might expect 
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Rayleigh-Taylor type instabilities for Cases 2,3, and 4, there 
is no body force which would drive the density difference; 
hence, the system will remain stable. 

This zero mean gravity state will be taken as the base 
state for the remainder of the investigations of this thesis. 


22 






Dispersion Solution 



(s/mo) poods uouB^Bdoid 
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Dispersion Solution 

air/silicone oil/air 
h=lcm, k=l cm ' 1 
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solid line: real root (stable solution) 
dashed line: imaginary root (positive values yield unstable solutions) 
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Effect of Wave Number on Stability 
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Effect of Wave Number on Stability 
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Effect of Height on Stability 
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Effect of Height on Stability 
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legend: height of middle slab, h (cm) 


CHAPTER 3 


MULTI-LAYER FLUID CONFIGURATION STABILITY IN THE PRESENCE 
OF A TIME-DEPENDENT PERIODIC ACCELERATION FIELD 

3.1 Problem description 

The results from Chapter 2 illustrated that for constant 

zero gravity the fluid configuration is stable for all 

parameter values. This, however, is a highly idealized case. 

A recent summary 22 indicates that the environment of board 

space shuttle is subjected to residual accelerations ranging 
-5 -3 

from 10 to 10 *g««rth at a frequency range up to 10 Hz. 

This chapter investigates the effect of periodic 

accelerations on the interface stability of a multi-layer fluid 

configuration without rigid boundaries. The configuration 

consists of a layer with finite height situated between two 

semi-infinite layers. All three layers extend to infinity in 

the horizontal plane. (See Figure 3.1). The accelerations 

are periodic about a zero mean gravity level, and are oriented 

normal (in the $ direction) to the interfaces. 

z 

The three fluids are inviscid, incompressible, 
irrotational, and immiscible. Surface tension is a property of 
the interfaces. The spatial dependence of the perturbation is 
considered to be wavelike. 

The dynamics of the above system are again governed by the 
continuity and Euler equations. These are non-dimensional ized 
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Configuration Geometry 


region 2 P2> $2 



region 1 Pb 0 I 


H 


z = F(t)e'( u+m y) 



Yiii 


region 3 P3» 4*3 

p = density of subscripted region 

= potential function of subscripted region 
y = surface tension of subscripted interface 


Figure 3.1 
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and linearized, resulting in a system of linear equations with 
time-dependent coefficients. 

A Floquet analysis is applied. The resulting system can 
be viewed as an eigensystem in the Floquet exponent. It is the 
value of this exponent which will determine the stability of 
the configuration. The linear stability of a perturbation to 
the interface is thus dependent upon six non-dimensional 
parameters: the density ratios of the outer to middle slabs, 
the Froude type number, the Bond type number at each interface, 
and the wave number. 


3.2 Equation development 

3 . 2a Governing equations and non-dimensionalization 
The governing equations are the continuity and the 
conservation of momentum equations for an incompressible 
and inviscid fluid. 


V*u - 0 


(3.1) 


au 

p + pu*7u 

at 


-7p - pG Q g(u f t)e 2 


(3.2) 
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where the time-dependence is apparent in the body force 
term. G q represents the peak value of the acceleration due to 
the periodic g-jitter. Equation (3.2) is to be linearized 
about a state of zero mean motion. Quadratically small terms 
are neglected (after expansion in a small disturbance 
parameter, c) . 

Expansions of pressure and velocity fields are as follows: 

p " p mean + cp '-* + *- (h.o.t.) (3.3) 

U ■ 0 + CU'..+.. (h.o.t.) (3.4) 

(Note the analysis considers zero mean motion (u -0).) 

'-mean 

The governing equations can be non-dimensionalized to 
yield: 

V-U' - 0 (3.5) 


P 5u' 

1 - -V p' (3.6a) 

P D ^ 


0 ■ -7 p 

*mean 



A 

e. 


(3.6b) 
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where u 


(3.7) 

(3.8) 

(3.9) 


= H(J f u 
x = H x 
t = (l/u f ) t 

P - P D H 2 Uj p (3.10) 

(Note that p D is the average of the density differences 
across each interface, and is the forcing frequency.) 

Equation (3.6b) is the mean conservation of momentum 

equation, and (3.5,3.6a) represents the perturbed system. 

Note that due to the periodic time-dependence, the mean 

pressure field will also be periodic in time. The parameter 
2 

(G q /Hu)j ) in equation (3.6b) is taken to be roughly of order 
one; this ensures the mean pressure to be of the same order. 

For convenience, the tildes will be omitted from this 
point forth in Chapter 3. All quantities are henceforth 
non-dimensional. Also, the primes will be dropped from 
perturbation quantities. 

The analysis is incompressible, inviscid, irrotational , 
and linear. A potential function can be utilized, giving rise 
to Laplace's equation. Perturbations are taken to be wave-like 
in the (xy) plane. The resulting differential equations (in z) 
are solved by separation of variables to yield: 
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^ = [Attle 112 + B( t) e -kz J e i(lxtm *> (3.11) 

0 2 = [c(t).- kz ] e i ( ix+ioy ) (J.12) 

* 3 - [d ( t) ® kz ] « i<lx+my) (3.13) 

A, B, C, and D are time-dependent coefficients, and k 
represents the wave number (k = 1 + m ) . ^ is the potential 

function. Subscripts indicate the region of interest. The 
pressure (both mean and perturbation) can be obtained from 
equations (3.6a) and (3.6b). 

A normal mode perturbation approach in the spatially 
independent variables is utilized; thus, the equations of each 
equilibrium interface can be written as: 

F.„- 2 - 1 - eE(t). 1 < 1X+B *> (3.14) 

F. m - 2 - 0 - cF(t) e*- (l x+m y( (3.15) 


where Fe defines the equation of each equilibrium interface, 
and E and F are time-dependent coefficients. 
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3 . 2b Boundary condition*; 


Three 

boundary 

conditions are 

imposed on each of 

the 

interfaces 

: (1) 

the kinematic 

boundary condition. 

( 2 ) 

continuity 

of the 

normal component 

of velocity, and (3) 

the 


normal force balance across the interface. 

The kinematic boundary condition states that at each 
interface 


- 0 
Dt 


(3.16) 


which can be expressed as: 


-=-= + u* V (Fe) - 0 (3.17) 

at 


Note that after linearization, the gradient of the interface 
equation has a contribution in the e -direction only. 
Imposition of this kinematic condition on each of the two 
interfaces yields the following: 


on Fe ZI , 


E (t) + k C(t)e“* - 0 


(3.18) 


on Fe Ij;i , -F (t) + k D(t) » 0 


(3.19) 
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A second boundary condition can be imposed in which the 
normal component of the velocity is continuous across the 
interface . 


a* d<t> 

on Fe TT , « (3.20) 

3Z 3z 


3 * 3 * 

on Fe , - (3.21) 

dZ 3z 


After a Taylor series expansion at each interface, the 
following relationships are obtained: 


A(t)e 2k - B(t) - -C(t) 


(3.22) 


A(t) - B(t) - D(t) (3.23) 

And finally, a (linearized) normal force balance across 
the interface is implemented. 


Pi - p » r 

^lower ^upper 


(3.24) 
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or, in non-dimensional form: 


„3 2 
H 


( p i 


ower ^upper 


= 7 


A 

n 


( 3 . 25 ) 


where r =* surface tension 

ft » linearized outward pointing normal to the 
interface 

Recall that the upper and lower pressures each have both a 
mean and a perturbation component. Contributions to the 
pressure at 0(c) (needed in equation (3.25)) involve both the 
perturbation pressure given in equation (3.26a) and a second 
term due to the wave itself. This second term is listed in 
equation (3.26b). 



(3.26a) 



cE(t) e 


second term 


i(lx+my) 


(3.26b) 



Note that the second term also contributes at the lower 
interface, with E(t) replaced by F(t). 

(Recall that primes have been dropped from the perturbation 
pressure . ) 

Substitution of equations (3.26a,b) into equation (3.25) 
for each respective interface results in the following 
relationships at 0(c): 


on Fe , 
II 



<P 2 -Pl> 

P D 


Fr*g(t) E(t) 


-MA(t)e K +B(t) 

k d 


e 


k ) 


+ 




- k 2 E(t) 


(3.27) 


on Fe 


III 


/ 



( Pj- p 3 ) 

P D 


Fr*g (t) F(t) 


+ 


— ( A ( t ) +B(t) ) 



= k 2 F(t) 


(3.28) 


where 


Bo 


p D H G o 


2,3 


Fr = 


II, III 


Hu). 


BO 2 j a nd Fr are Bond and Froude number type parameters. 
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, 3 . 23 ) , 


the C and D 


( 3 . 32 ) 



where 


P 


21 



31 



t 



Thus, the system has been reduced to four ordinary 
differential equations (in time) with four time-dependent 
unknowns: A, B, E, F. The time-varying forcing function is chosen 
to be periodic, with 


g(t) - cos(t) - | (e xt + e -lt ) 


(3.33) 


3.3 Application of Floquet theory 

Floquet theory can then be applied to system 

(3.29-3.32). This is done by expressing the four 

time-dependent coefficients as 

[A(t),B(t),E(t),F(t)] - n l!„tVB n ,E n ,F n ) « lnt < 3 ' 34 > 

where \ is the Floquet exponent. 

By substitution of equation (3.34) into the system (3.29- 
3.32), the four ordinary differential equations in time can be 
expressed as an infinite algebraic system with the Floquet 
exponent appearing as a parameter. 


44 



(X+in) 


( 1 + P 2 l )e 

P D1 


A n + (X+in) 


(P 21 -l)e 


-k 


D1 


B. 


n 


( 3 . 35 ) 


(P 21 -l) Fr 


D1 


< E n-l + E n+1> - 


Fr 


Bo. 


E n " 0 


(X+in) 


P D1 


+ (X+in) 


<^ p 3l) 

P D1 


B. 


n 


d-*31> 


Fr 


01 


< F n-l +F n+l> 



F 


n 


(3.36) 

0 


(X+in) E + kB e~ k - kA e k - 0 
n n n 


(3.37) 


(X+in) F„ + kB„ - kA„ - 0 
n n n 


(3.38) 


where n varies from -« to +®. This results in an infinite sys- 
tem of equations. The set of homogeneous equations 
(3.35-3.38) can be written in the form 
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(£-*§> X - 0 


( 3 . 39 ) 


which is the generalized eigenvalue problem. The Floquet 
exponent, A , acts as the eigenvalue. X is an infinite column 
vector containing the following terms: 



ib ~ AB) is the coefficient matrix of X. It is generated in 
groups of four rows corresponding to a particular value of n. 
Matrices £ and | are shown in Figures 3.2 and 3.3. 
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p21 - 1 



o 


o 



The solution to the single interface problem (Jacgmin and 
Duval) 12 , utilized Floquet theory and truncated at n = ’241. 
It was decided that the multi-layer system should be truncated 
at least at this level. In our analysis, n is truncated at 
1 25 1 . This generates a system of 204 equations in which the 
Floquet exponent is determined by eigenvalue methods. (A 
generalized eigenvalue routine is used) . Details of the 
algorithms are found in Kaufman (1974) 14,15 and Moler and 
Stewart (1973) 19 . 


3.4 Solution methodology 

To solve the large, sparse generalized eigensystem, a 
routine DGVLCG from the IMSL library package is utilized. 
(See Appendix 5.) This routine is based on the LZ algorithm 
described by Kaufman (1974), which in turn is similar to the QZ 
algorithm (Moler and Stewart , 1973) except that it uses 
elementary transformations whereas the latter uses orthogonal 
transformations . 

3.4a Preliminary checks 

Four checks were performed to verify the accuracy of the 


solutions. 



1) The system was converted to a standard eigensystem 

of form, (g - X J) X =* 0 , which was analyzed using DEVLCG 

from IMSL. This particular routine converts the matrix into 
a complex upper Hessenberg matrix, in which the eigenvalues 
are generated via the QR algorithm (Smith, 1976) 24 . There is 
agreement of the solution values. (See Appendix 6.) 

2) The original matrix was truncated at n * | 50 | , 

producing a much larger matrix. The same eigenvalues were 

obtained, with greater multiplicity of each root. 

3) The generated Floquet exponents were resubstituted 

into the linear system to compute the determinant by 
Gauss' method. The checks show our eigenvalues to 

be accurate. (See Appendix 7.) 

4) A limit case of the two interface system was 

investigated. In this case, P 21 “ 1.0 and Bo 2 » ». The physical 
interpretation of such a system is that the top and middle 
slabs have the same density, and their interface has zero 
surface tension. Hence, the system can be considered as a one 
interface configuration at Fe XII . 

To compare the results of the limit case, an analysis was 
performed for a one interface system following the methodology 
of sections 3. 2-3. 3. This new system reduces to two linear 
differential equations which are solved via Floquet analysis in 
the same manner as was done for the two interface system. 

(See Appendices 2 and 8.) A parameter, Q, appears with 
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Q = 1 . 0 * Bo 3 

(The constant, of unity exists due to fact that the Froude 
number identically equals one according to a comparison of the 
non-dimensional izations . ) 

A comparison of the two interface limit case with the one 
interface system for the same physical configuration is shown 
in Figures 3.4 and 3.5. The correlation between the two 
systems is evident. 

The results of these four separate checks provide great 
confidence in the numerical results. 

3 . 4b Solution interpretation 

The stability of the system can be determined by the 
numerical value of the complex eigenvalues. It is readily seen 
from equation (3.34) that the time-dependent coefficients will 
grow exponentially for positive real components of the Floquet 
exponents. Such a case will imply an instability of the fluid 
system. Thus, as the eigenvalues are generated, the presence 
of a single positive real part of A will dominate the system, 
causing it to be unstable. As this is a linear analysis, no 
information can be obtained concerning the finite amplitude 
(nonlinear) form of the configuration. 

The system of algebraic equations is non-dimensional. 
Thus, the linear stability of the fluid layers depends on six 
non- dimensional parameters: the two density ratios, the Froude 
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number, the Bond type number at each interface, and the wave 
number . 

A parametric study is performed to investigate the effects 
of parameter variation on the stability of the fluid system. 
The parameter space is defined by appropriate values of the six 
non-dimensional parameters pertinent to a microgravity 
environment . 

As stated previously, the complex Floquet exponent (A) is 
the eigenvalue of the system. The presence of a single, 

positive real component implies exponential growth of the 
interface and hence, a subsequent instability. Thus, we are 

concerned solely with the largest real component of the Floquet 
exponent. The values of this quantity are charted throughout 
the parameter space. A positive value of the largest real 
component of the Floquet exponent indicates an instability? a 
zero or negative value indicates stability. 
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3 . 5 Results 


To display most effectively the regions of instability, 
the largest (fastest growing) real component of the Floquet 
exponent is plotted as a function of the wave number for fixed 
values of Bond and Froude numbers and density ratios. 

The range for the parameters is as follows: 

0 . 1 < k < 5.0 

30 ^ = 1 . 0 , 0 . 1 , 0 . 01 , 0.001 

Bo 3 = 1 . 0 (Bo 2 ) , 2.0(Bo 2 ) 

Fr = 5.0, 1.0, 0.5, 0.1 

(Note that the quantities are nondimensional . ) 


These values correspond to physically realistic configura 
which might be expected In e mlcrcgravity environment as well 

as satisfying conditions for linear analysis. 

For each case, if the forcing function, g(t), were set to 
aero instead of cos(t), the configuration would be stable for 
all parameter space. The interfaces would simply oscillate 
with no growth of the amplitude. It is only with the forcing, 
and in the indicated parameter regions, that instabilities may 


occur . 

in Figures 3.6-3.25, the effect of Bo 2 <Bo 3 ) 
for different values of Fr and density ratios is 


on stability 
illustrated. 


_. - g. n -1 0 P =0.001225, Fr=5.0, 

For Fig. 3. 6, P 21 -1-u, P 31 

equal to Bo.,. The unstable wave number region 


and Bo 2 is set 
is broadest for 
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the largest Bo 2 (Bo 3 ) values. As Bo 2 (Bc> 3 ) is decreased, the 
unstable wave number region shrinks to encompass a smaller k 
range and tends towards lower k values. For low Bo 2 (Bo 3 ) 
values, the configuration is unstable to longer wavelength 
perturbations in the presence of periodic forcing. 

In Figure 3.7, Fr is reduced to 1.0 while other parameter 
values remain the same. Though the general qualitative trends 
follow, it is seen that the range of unstable wave numbers 

is broader than in Figure 3.6. 

Figure 3.8 shows an order of magnitude drop in Fr. Again 
as Bo 2 (Bo 3 ) is increased, the range of unstable wave numbers 
broadens. Likewise, as Fr decreases the unstable region 
encompasses a broader range of wave numbers. 

Similar results are elucidated in Figures 3.9—3.20, 
keeping parameter ranges the same for various density ratios. 
As Bo 2 (Bo 3 ) values are increased, the range of unstable wave 
numbers widens (corresponding to smaller wavelength 
disturbances) . 

The density ratios in Figures 3.21-3.23 pertain to a 
gas/liquid/gas configuration. The qualitative trends of 

varying Fr and Bo 2 (Bo 3 ) continue. However, the behavior is 
observed to be more peaky, with occasional regions of stability 
punctuating unstable wave number bands. For the lowest Fr 
(Figure 3.23), the unstable band shifts away from the low k 
region. Note that in Figure 3.21, the configuration is 
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stable when Bo 2 (Bo 3 ) is 0.001. Positive values of the real 
component of \ occur only for larger Bo 2 (Bo 3 ) values. 

Figures 3.24 and 3.25 show the effect of Bo 2 not 
equal to Bo^. In Figure 3.24 the Bo values are equal, whereas 
in Figure 3.25, Bo^ is twice Bo 2 , keeping all other parameters 
the same. Physically, the increase of Bo^ while keeping Fr 
fixed can be interpreted as decrease in the surface tension 
value at the lower interface. The dominant effect is to 
broaden the range of unstable wave numbers for each set of Bo 
values. Note that Figure 3.25 shows a narrow band of wave 
number stability near k=1.95. In general, the numerical value 
of the real part of the Floquet exponent (X) is increased for 
Bo^ twice the value of Bo 2 , indicating a faster growing 
"fastest growing" disturbance. 

The effect of varying Fr while holding Bo 2 (Bo 3 ) values 
fixed is illustrated in Figures 3.26-3.29. As Fr is decreased, 
the range of unstable wave numbers increases. Physically, this 
can be interpreted as a decrease in configuration stability for 
larger frequencies of the g-jitter. The behavior is typical 
for the various density ratios and Bo values which were 
considered. 

The effect of density ratio difference on stability is 
presented in Figures 3.30-3.34. Values of p 21 and p 31 
represent the density ratios of the upper and lower regions to 
that of the middle layer, respectively. Among cases indicated, 
the largest value of p Q1 [ = ( I p -l 1 + |p 3 -II ) /2 ] corresponds 
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to the case having the largest range of unstable wave numbers. 
Furthermore, it also corresponds to the largest values of the 
real component of the Floguet exponent. In general, as P D1 is 
decreased, the band of unstable wave numbers becomes more 
narrow. Results are typical and illustrative for the parameter 
space of concern. Note that Figure 3.34 compares 

three configurations of gas/ 1 iquid/gas with different gas 
densities. The three cases have similar results, indicating 
that the density of the gas layers is not too significant. 

In addition, the case in which both density ratios were 
set to unity, indicating .equal densities in all three regions, 
was addressed. In such a case, p Q1 equals zero; hence, 
Bo 2 (Bo 3 ) values are identically zero. Under the action of 
g-jitter, lack of density differences between the layers 
results in a stable configuration. One would expect the 
interfaces to merely oscillate in time. As a further check, 
the system was derived using a different definition of P D 
t p D=( p i+ p 2 +p 3 )/3 ]. Hence, Bo values are non-zero for equal 

densities in each region. Numerical results showed stability 
for all parameter values. 

The height of the finite middle slab is a physical 
quantity which appears in both the Bo and Fr nondimensional 
parameters. In particular, Fr is inversely proportional to the 
height, while Bo depends on the square of the height. An 
increase in height, H, implies a decrease in Fr and an increase 
in Bo. As was seen in Figures 3.6-3.25, an increase in Bo 
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corresponds to a larger region of unstable wave number space. 
From Figures 3.26-3.29, it is seen that the broadest region of 
unstable wave numbers occurs at smaller values of Fr. Thus, it 
is expected that an increase in H will result in a more 
unstable fluid configuration. This is confirmed in Figure 
3.35. Results are presented graphically for the case G q = 

_ A 

(10 *g«arth) , u f =0.1 Hz, and *ii =r m~ 50 dynes/cm. In 

addition, p 21 =0.8, and p 31 =1.2. Larger values of H 

correspond to a larger range of unstable wave numbers. 

The wave number at which the subharmonic (A=l/2) occurs, 
is plotted in Figures 3.36 and 3.37 for a range of Fr with 
given values of Bo 2 (Bo 3 ). It is seen that there is a gradual 
shift of the subharmonic to lower wave numbers (or longer 
wavelengths) as Fr is increased (ie., as the forcing frequency 
decreases). Figure 3.37 represents the case of unequal Bo 
values (Bo 3 =2*Bo 2 ) . Physically, this implies the surface 
tension of the lower interface is halved. It is seen that this 
case has subharmonics occurring at larger values of the wave 
number than the case of Bo 3 =Bo 2 (as shown in Figure 3.36). 

Stability boundaries of Fr versus k are plotted in Figures 
3.38 and 3.39. This is done for configurations of different 
density ratios (as indicated by the different area fill 
patterns) . Moreover, on each graph, multiple values of Bo 2 3 
are represented. The unstable regions are indicated by the 
rectangular "filled" regions. No meaning is ascribed to the 
width of the rectangles. In general, it is seen that an 
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increase in Fr while holding other parameter values fixed 
(corresponding to low frequency forcing) results in a smaller 
range of unstable wave numbers. Likewise, an increase in Bo 
values (corresponding to a decrease in the surface tension) 
relates to a broader region of instability, in terms of k. 

In general, fluid systems involving larger values of 
have wider bands of unstable wave numbers. This is evident in 
both Figures 3.38 and 3.39. Note also the "gaps" in the band 
of unstable wave numbers. These represent regions of stability 
of the fluid configuration. It is generally at higher k values 
(k > 1) that these bands of stability occur. 

The limit case (P 21 =1.0, Bo 2 =®) was used to compare the 
stability of the one interface configuration with that of 
the multi-layer fluid system. To compare the two, the 
parameters in each case are set equal at interface Fe^. In this 
way, the parameters are consistent in both problems. From 
Figures 3.40-3.41, it is readily seen that the multi-layered 
configuration is more unstable than the one interface fluid 
system. The range of unstable wave numbers is broader for the 
multi-layered case. In particular, note the contrast in the 
low k region (corresponding to large wavelengths) . That is, in 
the two interface system the very low k region is generally an 
unstable region as compared to the one interface case. In the 
one interface model, bands of instability are more frequently 
punctuated by narrow regions of wave number stability. These 
results are typical and illustrative for various density 
ratios . 


58 



Interface Model 



O <D 

c- <-> 

o 

00 Jr? 

•r-| QJ 

c 3 C 
D- 1— i 

E <N 

O 

U .ts 
£ 





V 

*2 m 
2 <N 
E <N 




mauodxa i3nbo[j jo msuodtuoD [esj 


59 


Q=1.0 Bo 2 =*>, Bo 3 =1 .0, Fr=1.0 



Comparison of 1 Interface Model 
with 2 Interface Limit Approximation 
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P2i=0 5 p2i=1.0, pu=0.5 

Q=0-I Bo 2 =~, Bo 3 =0. 1 , Fr=l .0 



Effect of Bo on Stability 
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Figure 3.6 
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legend: values of B02 (=803) 
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Effect of Bo on Stability 
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legend: values of B02 (-B03) 


Effect of Bo on Stability 
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legend: values of B02 (=803) 


Effect of Bo on Stability 



legend: values of B02 (= B03) 
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Figure 3.13 
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Figure 3. 14 
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Effect of Bo on Stability 
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Figure 3. 19 
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Figure 3.20 
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Figure 3.23 
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legend: values of B02 (=603) 


Effect of Unequal Bo on Stability 
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Figure 3.25 
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Effect of Fr on Stability 
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legend: values of Fr 
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Effect of poi on Stability 
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Figure 3.31 
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Figure 3.32 
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Effect of Height Variation 
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Subharmonic Case 
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Figure 3.36 


Subharmoni 
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Figure 3.3 


Stability Boundaries 
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Stability Boundaries 

lor given density ratios 
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Figure 3.39 
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Comparison of 1 Interface N 
with 2 Interface System 
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CHAPTER 4 


RESPONSE OF MULTI-LAYER FLUID CONFIGURATION TO SHORT-DURATION 
NON-PERIODIC TIME-DEPENDENT FORCING 

4.1 Problem description and stability considerations 

The residual accelerations which occur in a microgravity 

environment are generally time-dependent in nature. The 

special case of periodic g- jitter was addressed in Chapter 3. 

In addition to periodic forcing, residual accelerations may be 

of impulse type, due to such causes as station-keeping 

maneuvers and astronaut motion. This forcing, though 

non-periodic, would certainly still be time-dependent. 

This chapter investigates the effect of time-dependent, 

non-periodic accelerations on the interface stability of an 

idealized fluid configuration. The geometry of the system is 

the same as in Chapter 3 (see Figure 3.1). The accelerations 

are again oriented normal (in the e direction) to the 

z 

interfaces and have a zero initial value condition. 
Moderate-duration responses will primarily be investigated. 

three fluids are assumed inviscid, incompressible, 
irrotat ional , and immiscible. Surface tension is a property of 
the interfaces and is taken to be constant. The spatial 
dependence of the perturbation is considered to be wavelike. 
The fluid system is governed by the continuity and Euler 
equations. These equations are non-dimensionalized and 
linearized, resulting in a system of linear equations with 
time-dependent coefficients. 
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will be considered: 


1) An 


Two forcing functions 
exponential ramp followed by exponential decay, and 2) a ramped 
step function with both positive and negative values. These 
idealized functions were chosen to represent general impulses. 

An analytical approach is used to ascertain the asymptotic 
(mathematical) stability of the non-autonomous system. This 
analysis is presented in Section 4.3, with additional reference 
to Appendix 3 . Note that the system is non-autonomous due to 

the explicit appearance of time in g(t) . 

The system of first-order differential equations is 
integrated numerically utilizing Gear's stiff method 4 in order 
to solve for the time-dependent coefficients, which describe 
the response of the interfaces. The interface responses are 

plotted as a function of time. 

The presence of the "short duration" non-periodic body 

force functions do not modify the asymptotic stability of the 
multi-layer fluid configuration. Interest is in the level of 
system disturbance in the presence of the acceleration. Since 
the effects of viscosity are not incorporated into the 
analysis, the long time behavior of the system, predicted by 
the model, is not physically meaningful. 
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4.2 Equation development 


4.2a Governing equations and non-dimensionalization 
The governing equations remain those of conservation of 
mass and momentum for an incompressible fluid. Again, the 
analysis is inviscid. Linearization is performed about a state 
of zero mean motion. (Quadratically small terms are neglected 
after expansion in c.) The forcing function is of the form 

G 0 g(t) 

where g(t) is non-periodic, and may be represented by a ramp 
function, for example. G q is taken to be the peak magnitude of 
the forcing function. A positive forcing value is oriented in 
the negative e 2 direction. 

For this non-periodic forcing case, the non- 
dimensionalizations used are: 


u * vHgT u 
— o — 

X =■ H X 


t - / H/G t 
o 

p = P D HG o p 


(4.1) 

(4.2) 

(4.3) 

(4.4) 
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The non-dimensionalized disturbance governing equations assume 
the following form: 


7*u' = 0 


(4.5) 


au' 

— - -7 p' (4.6a) 

P D at 


0 


-7 p 


mean 


g(t) e. 


(4.6b) 


where tildes indicate non-dimensional quantities and primes 
denote perturbation values. Equation (4.6b) says the mean 
pressure instantaneously adopts a hydrostatic distribution with 
a magnitude governed by the instantaneous value of g(t) . 
Henceforth, the tildes will be dropped, and all quantities are 
to be considered non-dimensional. Also, the primes are 
dropped. 

As in Chapter 3, the pressure and velocity fields were 
expanded into a mean and perturbation component. The mean 
velocity, as stated previously, equals zero. 

An inviscid, irrotational , and incompressible analysis 
gives rise to a potential function ( u = 7 <p) , which can be 
substituted into equation (4.5) to yield Laplace's equation. 
This equation is solved in each region, yielding the same 
potential functions as in Chapter 3. 
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♦ x - [A(t)« kZ * B ( t ) e -kz 1 


! 1 e i(1>w 


(4.7) 


♦ 2 - fc (t) e~ kz ] e i(1 * +m >'> 


(4.8) 


* 3 = [D(t)e kz l e 1 ( ix+my ) 


where A,B,C,D are time-dependent coefficients. k represents 

2 2 2 

the wave number, where k * 1 + m . 

The interface shapes are identical to Chapter 3. Repeating 
them here, they are 


Fe XI = z - 1 - cE(t)e l(lx+my) 
Fe XII = z - 0 - cF(t) e 1 (lx+my) 


4 . 2b Boundary conditions 

Three boundary conditions are imposed at each interface: 
(1) the kinematic condition, (2) continuity of the normal 
component of the velocity, and (3) a normal force balance 
across each interface. 
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Imposition of the kinematic condition yields the following 


relationships between the time-dependent coefficients: 


on Fejj, 


-E - kCe” k * 0 


(4.12) 


on Fe IIIf -F + kD = 0 (4.13) 

The normal component of the velocity must match across each 
interface. Applying this condition gives the following 
relationships : 


on Fe XI , 


Ae 


2k 


- B = -C 


(4.14) 


on Fe___, A - B * D 


(4.15) 


The linearized normal force balance across each interface 
states : 


(h 2 g o ) 


'lower ^upper 


= 7 


A 

n 


(4.16) 


Note the upper and lower pressures each have a mean and 
perturbation component. Substitution into equation (4.16) 
follows the same methodology as that discussed in Chapter 3, 
resulting in the following system of equations: 
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B ( t ) 


-( 1 + P 2 l )e 

P D1 


A ( t ) 


( p 21~ X ^ e 
P D1 


-k 


+ 



21 

P 


-i)g(t) 

D1 



(4.17) 


^- p 3l) 

P D1 


A(t) 


( X + p 31> 

P D1 


B(t) 


(l-p 31 )g(t) 

P D1 



F (t) 


0 


(4.18) 


E(t) = k^A(t)e k - B(t)e“ k ] ( 4 . 19 ) 

F(t) - k[A(t) - B (t) 1 (4.20) 

Note that E (t) and F(t) represent the displacement of the 
interface due to the perturbation. 

This system has thus been reduced to four linear, 
non-autonomous differential equations of the form: 
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X =( g + 

in) 

rt 

IX 

(4 

• 21) 

where X is 

the unknown 

vector function 

containing 

the 

time-dependent 

coef f icients . 

g is a constant 

matrix and 

E(t) 


is the time-dependent matrix. (See Figure 4.1.) 

4.3 Asymptotic stability 

According to Sanchez 23 , non-autonomous linear systems of' 
equations (in the form of equation (4.21)) are asymptotically 
stable if three conditions are satisfied: 

(1) the characteristic polynomial of g is stable, 

(2) the matrix g(t) is continuous on 0 * t < ® , 

(3) /" II |(t)ll dt < ® . 

The functions which are selected will be shown to satisfy 
conditions (2) and (3). Condition (1) is ascertained by 
checking the characteristic polynomial of matrix g. For a 
stable solution of equation (4.21), the four roots of the 
characteristic polynomial of g must have non-positive real 
components (see Appendix 3) . Hence, solutions of equation 
(4.21) are bounded and stable. 

Having ascertained the asymptotic stability of the fluid 
configuration, the time response of the interfaces in the 
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5,, (l-p21 )g(0 6,i(p31l)g(0 
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- [(1 *p21 X» +f>3 1 )c* ♦ (1 p3i XI p2 1 )e 



presence of forcing is investigated. The "short duration" 
non-periodic body force functions which have been constructed 
do not modify the asymptotic stability of the multi-layer fl uid 
configuration. Interest is in the level of system disturbance 
in the presence of the acceleration, (eg. E(t) and F(t)). 
Since the effects of viscosity are not incorporated into the 
analysis, the long time behavior of the system is not 
physically meaningful. However, long-duration responses win 
be examined to investigate asymptotic stability. 

Note that for asymptotic stability, the forcing function 
must be bounded. That is, if g( t, was chosen to be periodic 
(say a cosine function), condition (3) of Sanches would be 
violated, and the system of equations would not be 
asymptotically stable. This does not imply that the fluid 
configuration is unstable to periodic forcing. In Chapter 3 , 
it was shown that there exist regions of parametric stability 
in the presence of periodic forcing. 


4.4 Results 

4 «4a Solution methodology 

The system of first-order differential equations, 
(4.17-4.20), is integrated via DIVPAG of the IMSL library. 
This routine utilises Gear's method to solve for the 
time-dependent coefficients*, (see Appendix 9 .) 
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It was found to be more convenient to represent the fourth 
order system in terms of <E,F,E,F). This is accomplished by 

differentiating equations ( 4 . 19 , 20 ) and substituting into 

equations (4.17,18) to eliminate A(t) and B(t) . This system is 
given as follows: 

6 - [ *W< ! 3> * *«4-*3>] F - l* lW< -, 7 , E (4.22, 


F " [ kfl 8 (S 2 +a 3 ) +k(fl 9* S 10 ) ] F ' Wj(« 6 -S 7 ) E 


where 

g {-U+P 31 )e k - e“ k (i-p 3l)) p oi 

(1+P 21^ 1+P 31 )ek + (p 21 _1) d“P 31 )e“ k 

-d+P 21 )g(t)e k (l+p 31 ) k 2 e k 

0 , - « _ 21 

(1-P 31 )Bo 3 


0 5 - g(t)e k 


D1 




V 2 n 

k P Dl e 
Bo 3 (l-p 31 ) 
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*6 



-i)g(t) 

D1 


*7 


k 


2 


Bo 


2 


0 


~2p 


8 


D1 


(l+p 21 ) (l+P 31 )e K + (p 21 -l) (l-p 31 )e 






k 2 p 


D1 


Bo 3 (l-p 31 ) 


^10 - 9(t) 


The system was integrated for specified non-zero (E q/ f ) 
values with E q/ F q both taken to be zero. The E(t) and F(t) 
coefficients define the time-development of the interface. e 

o 

and F o are determined from equations (4.12,4.13) where c =-0.05 

o 

and D q » 0.05. These values were chosen arbitrarily but are 
required to be small to satisfy the restriction to linearity. 
The signs vers chosen to ensure that both interfaces had the 
same direction for their respective velocity fields. This was 
decided solely to provide a more realistic physical system. 
Different values of C Q and D Q were investigated with 
qualitatively similar results. 
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Two time-dependent forcing cases (for t*0) are investigated 
and are defined as follows: 


if 0 s t < 1 
if 1 s t < 2 

if 2 s t < 3 

if 3 4 t < 5 

if 5 4 t < 6 

if 6 4 t < 7 

if t * 7 

exponential growth followed by 

2 represents a ramped step function 

with forcing in the positive and negative e directions. Each 

z 

forcing function is shown in Figures (4.11-4.34) as a solid 
line. Both cases have an initial value of zero forcing and are 
continuous on 0 4 t < • ; hence, condition (2) for asymptotic 
stability is satisfied. 

Condition 3) requires that the integral over infinity of 
the absolute value of the forcing function be finite. For 
forcing case 1 (exponential ramp) , the integral is solved by 
integration by parts. 


l: g(t) = te ( ~ t+1) 


2 : 


g(t) 


0 

t-i 

1 

-t+4 

-1 

t-7 

0 


Case l represents 
exponential decay. Case 


S 

0 


II te (_t+1) || dt 


< to 
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F or the step forcing function, case 2, the integral is 
subdivided into appropriate tine steps and is found to equal 4. 
(Note that due to the absolute value in the integral, negative 
values of the forcing do not cancel out positives.) Hence, for 
both forcing cases the integral is finite, and condition (3) of 
Sanchez is satisfied. 

The time interval extends in the range, 0 s t < t , where 
t is a value of time significantly greater than the outer 
bound of the forcing function. 

The time response of each forcing case is investigated for 
parameter values pertinent to a microgravity environment. 

k - 0.5, 1.0, 2.0 
Bo 2 - 1.0, 0.1 
Bo 3 - 1.0(Bo 2 ) , 2.0(Bo 2 ) 

These values correspond to physically realistic configurations 
where values of G Q range from lo" 3 to 10“ 5 * g**rth. 

Two fluid systems are considered: 

1) gas/ liquid/gas (P 21 »P 31 *0. 001225) 

2) liquid/liquid/liquid (p 21 *1.0, p 31 =1.5) 

The reverse of system 2), (p 21 »1.5, p 31 »1.0), is investigated 

since the forcing cases are directional, and the behavior is 
fundamentally different. 



4 . 4b Numerical results 

The shape of the interface is defined by E(t) and F(t) , 
wh ich are the perturbation displacements. The time response of 
E( t) and F ( t) are examined for the said parameter space as 
described in section (4.4a). 

Discussion of zero forcing: 

The case of zero forcing was studied to show the time 
behavior of E(t) and F(t) in the absence of transient 
accelerations. As expected, the numerical results were 
consistent with the fact that the configuration is 
asymptotically stable 23 . Results show an oscillatory pattern 
that neither grows nor decays exponentially in time (See 
Figures (4.2-4.10)). The perturbations are wavelike. Note 
that the variation of the interface perturbations in the zero 
forcing case is not uniform and sinusoidal. This is due to the 
coupling effect of the two interfaces, which have different 
velocities according to equations (4.12,4.13). with careful 
selection of constants C Q and D q for a specific wave number, 

the initial velocities of each interface could be set equal 

• • 

(E(t)-F(t)). in the subcase of E(t)-F(t) and p -p., , the 
time variation of the perturbations is sinusoidal. Although 
equal interface velocities provides a more uniform wave on the 
interfaces, it is recognized as a special case. The general 
case, with fixed C Q and D Q and hence unequal E and F, is 
investigated in the results. 
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Also note that as the wave number increases (Figures 
4 . 2, 4. 3, 4. 4), the amplitude of the perturbation becomes larger 
in magnitude. This is due to the k dependence of the interface 

velocities (equations 4.12,4.13). Moreover, note that the 

\ 

selection of C Q and D q must be such that all quantities do not 
violate linear theory. 

For a fixed wave number for the zero forcing cases, the 
amplitudes of E and F are smaller for the gas/liquid/gas 
configuration than for the liquid/liquid/liquid systems 
(compare figure 4.9 with 4. 3, 4. 6 where k-1.0). This is due to 
decreased dynamical effects from the gas regions. Different 
wave numbers provide similar results. 

Discussion of impulse forcing (exponential ramp) : 

Figures 4.11-4.28 show moderate-duration responses of each 
interface to impulse forcing for the specified parameter space. 
Figures 4.11-4.19 correspond to the exponential ramp forcing 
case. Note that the forcing function is displayed on the 
graphs as a solid line. Physically, this function simulates 
short-duration impulse forcing which might be due to 
disturbances such as astronaut motion. By no way do the 
selected functions represent the entire class of possible 
impulses. It should be pointed out that positive values of the 
forcing correspond to accelerations in the negative e 

z 

direction. 


112 



Figure 4.11 represents a liquid/liquid/liquid configuration 
with the most dense fluid being on the bottom (p 21 =l.o, 
p 31 =1.5). The interfaces oscillate in time with a fairly 
periodic motion. It is clear that the response is greater in 
magnitude for higher values of Bo 2 (Bo 3 ). Bo values are 
inversely proportional to surface tension; hence, an increase 
in Bo is associated with a decrease in the restoring force at 
the interfaces. A decreased restoring force will lead to 
enhancement of the interface displacement. This trend is 
typical throughout the results. Note also that there is an 
increase in the period of the perturbation for higher Bo 
values. Henceforth, discussion will pertain to responses for 
Bo 2 (Bo 3 )=1.0 unless otherwise noted. 

Comparing Figure 4.11 with 4.2 (the same fluid system with 
no forcing, k*0.5), there is an enhancement of the interface 
displacements in the presence of forcing. The upper interface 
in the zero forcing case has AE -0.13 which increased to 

MX 

0.21 in the presence of the ramp forcing. This is a 62% 
amplification due to impulse forcing. Likewise, at the lower 
interface, AF increased from 0.14 to 0.24 for a 72% 

MX 

enhancement. (Note: AE ■ E - E ) . 

MX MX Bin 

The wave number in Figure 4.12 is increased to 1.0 while 
holding other parameters fixed. Oscillatory motion still 
occurs. Comparing with Figure 4.3 (zero forcing for the same 
configuration), there is no enhancement at either interface. 
In Figure 4.13 (k=2.0), the ramp forcing results in a smaller 
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interface displacement as compared to that of the zero forcing 
case (Figure 4.4). Recall that positive values of g(t) 
correspond to accelerations in the negative e ( ’’downward" ) 
direction. The greatest forcing for the ramp function occurs 
between t=0 and t=5. During this time period, the interfaces 
for Bo-1.0, particularly the lower interface (F) , are 

accelerating in the positive e direction. The net 

acceleration is smaller than for the zero forcing case, hence 
the "negative enhancements" of -30% and -40% at the upper and 
lower interfaces, respectively. 

Note that in both Figures 4.4 and 4.13, for Bo-0.1 and 
k— 2.0, both interfaces are accelerating in the negative e 

z 

direction during the time period in which the peak of the ramp 
forcing occurs. Note again that the ramp function is in the 

A 

negative e 2 direction. Thus the net acceleration is greater 
for the forced case with amplitude enhancements at upper and 
lower interfaces of 15% and 30%, respectively. 

Note that the results are contrary for Bo»1.0 and Bo=0.1 
at this k value of 2.0. 

The configuration is inverted for Figures 4.14-4.16, that 
i 9 • P 21 " 1 * 5 ' Such a configuration has an unfavorable 

density gradient with respect to the direction of the forcing 
(ie., a more dense fluid is oriented above a less dense one). 
The enhancement of the perturbation amplitude due to forcing 
(Figure 4.14) compared to the zero forcing case (Figure 4.5) 
£° r k-0.5 is 67% at the upper interface (E) and 90% at the 
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lower interface (F) . For k-l.o, the upper and lower interface 
amplifications due to ramp forcing are 12% and 11%, 
respectively. As with Figure 4.13 (p 21 =i.o, p 31 =i.5), the 
situation is reversed for k»2.0, with the non-forced case 
exhibiting larger perturbation amplitudes. Inspection of the 
interfaces for zero forcing (Figure 4.7) at Bo=1.0 show 
perturbation accelerations in the opposite direction of the 
impulse forcing during the critical time period (t=0 to t=5) . 

In general, the enhancements for the unfavorable density 
gradient configuration (p 21 «1.5, P 31 »1.0) are greater than 
those of the favorable one. 

Figures 4.17-4.19 correspond to a gas/liquid/gas 
configuration. Difference between forced and unforced cases 
are most dramatic at k-2.0. In Figure 4.19, there is "negative 
enhancement" for the ramp forcing as compared to the zero 
forcing (Figure 4.10). The upper interface exhibits little 
change, but the lower interface (F) has a 31% decrease in 
perturbation amplitude in the forced case. As in the 
liquid/liquid/liquid configurations, the cause is an interface 
acceleration in the opposite direction of the impulse forcing. 

In general, in the presence of ramp forcing the enhancement 
is greater at the lower interface. Recall that the forcing is 
mono-directional (downward) ; hence there is a true "upper" and 
a "lower" interface in terms of the acceleration field. 
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Discussion of impulse forcing (bi-directional step) : 

The step forcing function is employed in Figures 4.20-4.28. 
This implies that the acceleration is bi-directional (that is, 
an interface may experience a favorable density gradient at one 
moment and an unfavorable one at another time in its history) . 

Figures 4.20-4.22 correspond to a liquid/liquid/liquid 
configuration (Pj^-l.O, p 3 ^*1.5). Again it is observed that 
high surface tension (low Bo values) relate to minimal 
distortion of the interfaces. As with the exponential ramp 
forcing function, there is enhancement of perturbation 
amplitude at low wave numbers as compared to the zero forcing 
cases. In Figure 4.20 (k-0.5) the upper and lower interfaces 
are enhanced by 160% and 200%, respectively. This 
amplification due to the step function is considerably larger 
than the enhancement of the same configuration in the presence 
of ramp forcing (Figure 4.11). For k»1.0 (Figure 4.21) there 
is slight enhancement of 3% and 10% for interfaces E and F. 


Figure 4.22 

(k-2.0) 

shows "negative enhancement". Note 

on 

Figure 4.22, 

the lower interface response 

for 

Bo-1.0 is 

in 

"phase" with 

i the 

step forcing function 

but 

opposite 

in 

direction. 

During 

the period of forcing 

(t-1 

to t— 7 ) , 

the 

perturbation 

is in 

effect opposed, hence 

the 

reduction 

in 


amplitude. The long wavelength perturbations (low k values) 
are not M in phase" with the selected forcing functions; 
therefore, enhancement occurs even if the interface 
acceleration is opposite in direction to the forcing. 
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Figures 4.23-4.25 are for the case with the density ratios 
reversed and offer qualitatively similar results. since the 
configuration is subjected to bi-directional forcing, there is 
no true "upper" or "lower" interface with respect to the 
orientation of the forcing. Figure 4.23 (k-0.5) shows 

enhancements of 215* and 170* for interfaces E and F, 
respectively. For k-1.0 (Figure 4.24), the amplification of E 
as compared to the zero forcing case is 17* and for F, 21 *. 
For k-2.0 (Figure 4.25), there is "negative enhancement" at 
both interfaces. Again, note the "phasing" between the lower 
interface response and the period of the forcing. 

It should be pointed out that the greatest enhancement in 
perturbation amplitude for the bi-directional forcing occurs at 
the interface with P 21 (P 31 )-1.5 rather than l.o. That is, the 
interface with a density difference across it experiences 
greater enhancement. In general, this trend is typical for 
Figures 4.20-4.25. 

The gas/liquid/gas configuration in the presence of step 
function forcing is represented by Figures 4.26-4.28. m 

general, such configurations have more uniform oscillations. 
As occurred with exponential ramp forcing, it is the k=2 . 0 
configuration in which differences between the forced and 
unforced case are most dramatic. For k-2.0 in the forced case 
(Figure 4.28), there is an increase in the amplitude of the 
interfaces: E (25%) and F (32%). Recall that for all other 
cases with k-2.0, there is “negative enhancement". The 
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"phasing” between the interface response and the forcing period 
is the critical factor. Note on Figure 4.28, when the forcing 
is clipped (t-7), both interfaces are perturbed at a large 
amplitude, hence a "positive enhancement". Comparing with 
Figure 4.25, both perturbations have small displacements at 
t-7; the oscillation is in "phase" with the forcing period. 
Therefore, this case has "negative enhancement". At low k 
values (long wavelengths), for the selected forcing functions, 

phasing does not occur; thus in general, the impulse forcing 
enhances the perturbation. 

The effect of unequal Bo values is addressed in Figures 
4.29-31, where in each case So 3 -2*Bo 2 . That is, the surface 
tension of the bottom interface is half that of the top. Note 
that the case of equal Bo values appears on the graphs for 
purposes of comparison. The cases involving the doubling of 
one of the Bo values show an enhancement of the interface 
displacement as compared to the equal Bo case. in Figure 
4.29, the amplification for Bo 3 -2*bo 2 is 281 greater at E and 
100% larger at F as compared to the equal Bo case. Figure 4.30 
Shows a different configuration but similar results. The effect 
of doubling BOj actually has a slight "negative enhancement" 
(-1%) on Interface E, but a loo* amplification on F. Figure 
4.31 represent, a gas/liquid/gas configuration subjected to 
ramp forcing. Doubling Bo, corresponds to amplifications of 
26% on interface E and 58% on interface F. In all cases of 
Bo 3 -2*Bo 2 , the greatest enhancement occurs on interface F. 
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Doubling the Bo value on this interface corresponds to halving 
its surface tension. it is of interest that in the 
liquid/1 iquid/1 iquid configurations (Figures 4 . 29 , 4 . 30 ), 

doubling the Bo value of an interface corresponds to a loot 
amplification of that interface. 

It has been discussed that consideration of "long-duration" 
responses is not physically realistic due to the absence of 
viscosity. surely viscous effects would play an important 
damping role (ie., in time, we would expect that long after the 
forcing dies out, the perturbation amplitude should be damped) . 
For an inviscid system the perturbations, even in the absence 
of any forcing, continue to oscillate ad infinitum. Therefore, 
there is no physical relevance to the long-duration response. 
However, the extension of calculations out to these larger 
times yielded numerical results which are consistent with the 
known (mathematical) asymptotic stability of the configuration. 

Figures 4.32-4.34 display long-duration responses for 
various parameters. The two figure, in each left column 
represent zero forcing. The oscillatory nature of the 
perturbation is apparent. The two figures in the right column 
represent the same configuration in the presence of the 
designated forcing. Although all three figure, show an 
enhancement in the amplitude in the case of forcing, this is 
not to suggest that impulse forcing causes enhancement of the 
zero forcing case for all parameter space. (Recall the ramp 
forcing case for k- 2.0 where negative enhancements occur.) The 
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period of the oscillation in the forced cases is approximately 
equal to the period of the zero forcing case. During this time 
period (0<t<200) , the response of the interfaces is not growing 
exponentially in time. Nondimensional time periods up to 
t =1000 were examined with similar results. 

In general, the presence of impulse forcing causes 
enhancement of interface displacement (in the case of low k, 
long wavelength disturbances) . Depending on the phasing 
between the oscillation of the interface and the period of the 
forcing function, a reduction in interface amplitude may occur 
for some parameter space, particularly at higher wave numbers. 

The interface behavior for a given configuration may be 
very different for various impulse accelerations. Recall that 
the displacements for the gas/1 iquid/gas configuration at k=2.0 
were smaller in the presence of ramp forcing but were enhanced 
when subjected to step forcing. The possibility of enhancement 
could cause adverse effects on materials processing 
applications. 
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Time Response of the Interfaces 
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Figure 4.3 
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Time Response of the Interfaces 
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Time Response of the Interfaces 
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Time Response of the Interfaces 
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Time Response of the Interfaces 
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Figure 4.20 
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Time Response of the Interfaces 
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Figure 4.21 
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Time Response of the Interfaces 
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Figure 4.22 
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Time Response of the Interfaces 
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Time Response of the Interfaces 
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Time Response of the Interfaces 
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Time Response of the Interfaces 
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Time Response of the Interfaces 
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Time Response of the Interfaces 

effect of unequal Bo 






Long Duration Response: Comparison of Zero Forcing with Impulse Forcing 
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Figure 4.34 



CHAPTER 5 


CONCLUSIONS 

The effect of microgravity environment accelerations on 
the behavior of a multi-layered idealized fluid configuration 
has been investigated. The analysis was linear, and each fluid 
region was considered inviscid, incompressible, irrotat ional , 
and immiscible. A normal mode approach was taken with regard 
to the spatial variables. 

As a preliminary study, the stability of the configuration 
was investigated in the presence of constant acceleration 
fields. Dimensional equation development resulted in a 
dispersion relation. The nature of the roots of the dispersion 
relation determined the stability of the configuration. 

Three parameters were studied, including the wave number 
of the perturbation, height of the middle layer, and the value 
of the constant forcing. The stability regimes of these 
parameters were investigated for various configurations 
involving air, water, and silicone oil. 

The results show that the configuration is most stable to 
larger values of the wave number. This implies that the fluid 
system is susceptible to long wavelength perturbations. The 
change in height of the middle layer has a negligible effect if 
the quantity h/X is greater or on the order of one. For values 
of h/X < 0(1), faster growing instabilities are associated with 
smaller values of h/X. 
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value is 


Results indicate that as the constant forcing 
decreased, the configuration becomes more stable. The limit 
case of zero forcing was investigated and found to be stable 
for all parameter values. The zero mean gravity state served 
as the basis for the ensuing time-dependent cases. 

For the periodic case, the equations were 

non-dimensionalized. Floquet theory was applied to the system 
of equations (3.29-3.32) resulting in an infinite set of 
algebraic equations. A truncation was made, and the problem 
was posed as a generalized eigenvalue problem. Solutions to 
the eigensystem determined the stability of the configuration. 

Six non-dimensional parameters were investigated: the 

Bond type number at each interface (Bo 2 ,Bo 3 ), the density 
ratios of the outer to middle layer (02i'^3i^ ' th® Froude type 
number (Fr) , and the wave number (k) . Ranges of values studied 
are pertinent to a microgravity environment and satisfy 
conditions of linearity. Results indicate several trends 
involving these parameters. 

As Bo values are increased, the configuration becomes more 
unstable. That is, the unstable range encompasses a wider 
range of wave numbers. This trend can be interpreted as 

resulting from a decrease in the surface tension (inversely 

proportional to Bo) at the interfaces which expectedly would be 

more unstable. For unequal Bo values (Bo - 2*Bo.) , the range 
of unstable wave numbers is even greater. 

A density difference parameter (p^) is expressed in 
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terms of p 21 and p 31 - Essentially it equals the average of the 

density differences across each interface. in general, as o 

k di 

increases, the configuration becomes more unstable. This trend 
corresponds within certain "families" of configurations (for 
example, gas/liquid/liquid or liquid/liquid/liquid, etc.) 

An decrease in the Froude type parameter (inversely 
proportional to the square of the forcing frequency) 
corresponds to larger bands of unstable wave numbers. Hence, 
the configuration is more unstable to high frequency forcing. 

The multi-layered fluid system was found to be "more 
unstable" than the one interface configuration. That is, the 
range of unstable wave numbers is smaller in the one interface 
case. In particular, one area of contrast was in the very low 
k region where regions of stability were present for the one 
interface case. 

For the non-periodic forcing case, the non-dimensional ized 
equations resulted in four ordinary differential equations in 
time. The system was integrated numerically, and the time 
responses of the interfaces were obtained. 

Two short -duration impulse type functions were imposed on 
the system. Asymptotic stability of the fluid system in the 
presence of short-duration accelerations was ascertained via 
mathematical analysis, and the numerical results were 
consistent. The interfaces respond in a wavelike fashion, but 
do not grow exponentially in time, providing that certain 
conditions on the forcing function are satisfied. 
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In general, the presence of an impulse causes enhancement 

interface perturbation amplitudes (for long wavelength 
perturbations) as compared to the zero forcing case. For 
higher wave numbers, different impulse accelerations can affect 
a given configuration quite differently. if the oscillation of 
the perturbation has "phase correspondence" with the period of 
the forcing, a reduction in interface amplitude may occur. 
Perturbation enhancement is generally greater in the presence 
of the bi-directional step forcing as compared to the 
one-directional ramp forcing. While the wave is not growing 
exponentially in time, enhancement could cause undesired 
consequences for an experiment. For example, a solidification 
experiment could be adversely affected by the presence of 
impulse forcing. 

The results of the idealized fluid system are 
qualitatively relevant to specific configuration geometries. 
For example, it was determined that the multi-layered case is 
generally unstable for low wave numbers (long wavelengths) . 
Certain float zone processing techniques involve a fluid column 
which is multi-layered. Such a configuration would need to 
avoid long wavelength perturbations. In general, it was found 
that the multi-layered configuration has a wider band of 
unstable wavelengths than the single interface fluid system. 
Hence, any space-processing geometry involving multiple layers 
of fluids would be more susceptible to instabilities than a one 
interface configuration. 
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Additionally, the subharmonic 


case 


is relevant to 
space-processing applications. It was discovered that the 
fluid system is most unstable at low values of Fr (inversely 
proportional to the forcing frequency) . The investigation into 
the subharmonic case showed that at low Fr values, the 
subharmonic (A=l/2) occurs at higher wave numbers. 


This study involves values of non-dimensional parameters 
which are relevant to a microgravity environment. 
Configurations involving fluids of specific interest may be 
investigated. For example, a typical configuration may have 
the following dimensional parameters: p^= 0.8 g/cra 3 , 

*IU* 25 dynes/cm, G q » 10 3 * g e «rth, 0.5 s” 1 , and H= 4.0 

cm. These values, according to the definition of the non- 
dimensional parameters, correspond to values of Fr=l .54 and 
b °2 f b °3 ) *° • 02 • The configuration parameters are typical of 
what may be expected in microgravity processing applications. 
Fluid systems of specific interest may be investigated in such 


a manner. 

The multi-layer configuration utilized in this study was 
idealized. In an actual space-processing application, the 
fluid system would be bounded in space; the boundary 
conditions pertinent to the container would need to be 
considered. A suggested area of future investigation is to 
consider finite configurations. 
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appendix 1 


Utilization Of WVNET Rssnnrr^ 


The numerical solution and graphical representation of the 

present analytical problem retires a well-integrated host of 

computer resources. The CMS system was accessed via WVNET on a 

VT320 terminal. A remote site at the Engineering Sciences 
building was used. 


Al.i Numerical Results 

The numerical results for Chapters 2, 3, and 4 were obtained 
by accessing several routines from the IMSL library’? Programs 
which were utilized are found in Appendices 4-9. One solution, 
in the case of periodic forcing, involved the eigenvalues of a 
ry large complex matrix system. An enormous amount of 
storage space was required for computation. upon request 
WVNET increased the storage capacity from 4M to 12M. This was 
sufficient to run the programs. Alternatively, temporary disKs 
could be accessed to provide the necessary space. The 

following step, were taXen to declare the temporary disk space: 
TDSK 192 DISK B CYL 15 


FORMAT 192 B 
RELEASE A 
RELEASE B 
ACCESS 192 A 
ACCESS 191 B 
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space . 


The 


Thesa steps free 15 spare cylinders of disk 
computer now interprets this disk as the A disk and the 
original disk as the B disk. Hence, to bring files over to the 
temporary dlskspace, the following command must be used: 

COPY filename FORTRAN B filename FORTRAN A 


Once a data file is created, the file can be transferred back 
to the permanent storage using the following command: 

COPY filename filetype A filename filetype B 


This file is now saved in the permanent directory. After 
logging off, the temporary disk memory will be destroyed. This 
method was solely used prior to the increase of storage space. 

A typical session using the expanded memory is as follows: 
(After logging on to CMS via WVNET.) 

DEF STOR 12M 
I PL CMS 
GETDISK IMSL 
F0RTVS2 filename 

GLOBAL TXTLIB VSF2F0RT CMSLIB IMSL1 IMSL2 
GLOBAL LOADLIB VSF2L0AD 
LOAD filename 
START 
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Execution will create the desired datafile. 


A1.2 Graphical Results 

Two options were explored for graphing the results. 
Initially the data was downloaded to a diskette via Kermit, 
which in turn was plotted using Lotus/123 graphics package on a 
Zenith DS computer. While the output was satisfactory, it was 
inconvenient and time-consuming to change terminal sites. 

The second, and preferred, option was to access CMS 
directly through a WVNET line connected to a Macintosh II PC . 
This was accomplished via VersaTerm and VersaTerra Pro 
communications. The program calling IMSL routines was run in 
the same manner as with a VT320. The data was then transferred 
to a SAS/Graph routine emulating TEK4014 device, which 

presented the results graphically. A typical graphing session 
is as follows: 

COPY datafilename filetype A FOR017 LISTING A 

SAS filename of sas program 

TEK4014 


A typical SAS/Graph program is as follows: 

CMS FILEDEF FOR017 DISK FOR017 LISTING; 
DATA ; 

INFILE FOR017; 
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INPUT X Y BO; 

PROC GPLOT ; 

PLOT X * Y = BO ; 

SYMBOL1 I=SPLINE L=l; 

SYMBOL2 I=SPLINE L=21 ; 

SYMBOL3 I=SPLINE L-20; 

SYMBOL4 I=*SPLINE L-22; 

This routine will take three columns of data as input and 
graphically sort according to equal values of Bo. 

The plots are converted to MacDrav files from which 
hardcopies were obtained from MacDraw i and II graphic 
packages. The advantages to this option are the one-terminal 
site capabilities as well as good resolution. 
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Utilization of WVNET Resources: 


Method 1: 



A Zenith DS PC 


data downloaded to 
diskette via Kermit 


graphical results 
using Lotus/123 
graphics package 


- satisfactory output 

- inconvenient and time-consuming to change sites 
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method 2: 



- option of choice 

- one-terminal site capabilities 
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APPENDIX 2 


One Interface System of Equations 

The fluid system and analysis is the same as for the 
multi-layer configuration except that there is only one 
interface. The upper region is subscripted by a 2 and the 
lower region is subscripted by a 1. 

The same governing equations are utilized with the 
following non-dimensional izat ion : 


«5 0 /w f ) u 

(A2.1) 

2 


V w f ) x 

(A2 . 2) 

U/w f ) t 

(A2 . 3) 

W w f > 2 p 

(A2.4) 


The interface is given by: 

Fe - z - 0 - eC(t) e i ( lx+m y) 
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The velocity potentials in each 


region are described by: 


h - tA(t)e kz , e* (Ix+my) (A2 . 6) 

* 2 - C B ( t ) e" kz ] e i<l* + ny> (A2>7) 

(Note that tildes have been dropped. All quantities are 
non-dimensional . ) 

Application of boundary conditions is similar to that of 
the multi-layer configuration. This system reduces to two 
linear differential equations which are solved using Floquet 
analysis. The one interface system is as follows: 


(A+in) C - kA » 0 
n n 


(A2 . 8) 


(A+in)A n - 


2(1+P 21 ) 


(C n-l + W + 


k 2 p 


D1 


Qd+P 21 ) 


n 


(A2.9) 


where 


( 1 “ P 2 1 ) 


D1 
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3 


Q 


P 01 G O 





( A2 . 10) 


but there is no H, thus the lengthscale (G o /u f 2 ) is used. 


Q - Bo (1.0) 2 


(A2. 11) 


The problem can now be posed as an eigensystem which is 
truncated and solved numerically. 
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APPENDIX 3 


Stabili ty of Characteristic Equation 


According to condition (l) of Sanchez 23 (Section 4.3), the 
characteristic polynomial of § must be stable (ie. the four 
roots of the polynomial must have non-positive real components) . 
The polynomial is of the form 


A 4 + aA 2 + b 


or 


as + b = o 


(s=A 2 ) 


(A3 . 1) 


y 1 2 

-a±a-4b 

s » 

2 


(A3 . 2) 


For guaranteed non-positive real components: 

i) a a 0 

ii) b > 0 

iii) a 2 > 4 b 

Conditions i) and ii) are satisfied simply by the signs of 
their components. Condition iii) requires that the following 
inequality must hold true 
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(A3 . 3) 


T i ( r o + ^) 


T 4 T 5 T 5 


where 


. „ ~ k3 Ppi 

1 -(l+p 21 ) (l+p 31 )e k + (p 21 -l) (l-p 31 )e" k 


X 2 " 


-(l-p 31 )e“ k + (l+p 31 )e k 

B °2 P D 1 


d“P 21 )« * + (1 + P 01 )« 


21 ' 


*>3 *01 


«*D1 

4 -(1+P 21 ) (l+p 31 )e k +(P 21 -D (1-P 31 )e" k 


V «v 

t 5 - ( - e K + e K ) 


k 


6 


Bo 2 ®°3 
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As an analytical check, condition iii) was investigated 

for p ■ P 21 31 P31' B°2 = ®°3* The following requirement of 

stability was obtained: 

0 > - 4 p - (l-p) 2 e~ 2k ( A 3 . 4) 

This is true for all values of p. 

Using a root finder, the roots of the characteristic 
polynomial were determined for the other parameter cases, and 
all roots had non-positive real components. Hence, condition 
( 1 ) of Sanchez is satisified. 
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APPENDIX 4 


Dispersion Solution - Chap 2 

This program solves the dispersion relation which was 
derived in Chapter 2 of the thesis (equation (2.19)). The 
dispersion relation is a fourth order polynomial, the roots of 
which are the propagation speeds of the disturbance. 

An IMSL routine, 'DZPORC', is called. This routine is a 
complex root finder. Solutions are obtained for the various 
configurations of interest. 
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C THIS IS IN FILE: 'LAYERS FORTRAN' 

C 

C THIS ROUTINE SOLVES THE DISPERSION RELATION FOR THE CASE OF 
C CONSTANT GRAVITATIONAL FIELD. 

C THE DISPERSION RELATION IS A FOURTH ORDER POLYNOMIAL. THE 
C FOUR ROOTS ARE COMPLEX AND ARE SOLVED BY CALLING A ROOT 
C FINDER ROUTINE, 'DZPORC', FROM THE IMSL LIBRARY. 

C 

C SOLUTIONS ARE OBTAINED FOR THE FOUR CONFIGURATION CASES ACROSS 
C THE PARAMETER SPACE OF INTEREST. 

C 

C 

c 

IMPLICIT REAL* 8 (A-H,0-Z) 

DIMENS ION DEN ( 3 ) , GAM ( 3 ) 

REAL*8 COEFF (5) , ACFS ( 5) 

COMPLEX* 16 ROOT (4) 


COMMON/ DAT/ R 1 , R2 ,R3 , AH, AG 2 , AG3 , AGRAV, AWN 

NDEG-4 

GO-980. 0D0 


OPEN (UNIT- 14 , STATUS- ' NEW ' , FILE- ' FORO 14') 


DENSITIES OF FLUIDS 
DEN ( 1) -1 . 0D0 
DEN (2 ) -0 . 96D0 
DEN ( 3 ) -0 . 00122 5D0 

SURFACE TENSIONS 

GAM(l) — 72 . 0D0 
GAM (2) -40. 0D0 
GAM ( 3 ) — 2 5 . 0D0 


DO CASES 


CASE 1: AIR/SILICON OIL/WATER 
R2-DEN(3) 

Rl-DEN (2) 
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c 


c 

c 

c 

40 

c 

c 

80 

90 

100 

CASE 
C 


R3=DEN ( 1 ) 

AG2=GAM ( 3 ) 

AG3-GAM ( 2 ) 

WRITE ( 14, *) 'AIR OVER SILICON OIL OVER WATER' 
WRITE ( 14, *) ' DEN2- ' , R2 
WRITE (14,*) 'DEN1«',R1 
WRITE (14,*) 'DEN3-' ,R3 

WRITE (14,*) 'SUR TEN2=' , AG2 , 'SUR TEN3=',AG3 

DO 100 13-1,5 
DO 90 12-1,6 
DO 80 11-1,7 

AH-0 . 5D0+ (13-1) *0 . 25D0 
AGRAV-GO* (1.0D0-(I2-1) *0.2D0) 

AWN-0. 25D0+(I1-1) *0 . 5D0 

WRITE ( 14 , * ) 'H-' , AH, ' GRAV- ' , AGRAV 
CALL DISP(ACFS) 

DO 40 1-1,5 
COEFF ( I ) -ACFS ( I ) 

CONTINUE 

CALL DZPORC(NDEG, COEFF, ROOT) 

WRITE ( 14 , *) 'RT1-' , ROOT ( 1 ) 

WRITE (14 , *) ' RT2- ' , ROOT ( 2 ) 

WRITE ( 14, *) ' RT3- ' , ROOT ( 3 ) 

WRITE ( 14 , * ) ' RT4- ' , ROOT ( 4 ) 

CONTINUE 

CONTINUE 

CONTINUE 


: AIR/WATER/AIR 
R2-DEN(3) 

Rl-DEN(l) 

R3-DEN ( 3 ) 

AG2— GAM(l) 

AG3-GAM(1) 

WRITE (14 , *) ' ' 

WRITE (14, *) ' ' 

WRITE ( 14, *) 'AIR OVER WATER OVER AIR' 
WRITE ( 14 , * ) ' DEN2-' , R2 
WRITE ( 14 , * ) ' DENI- ' , R1 
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WRITE ( 14 , * ) ' DEN3- ' , R3 

WRITE (14,*) 'SUR TEN2= ' , AG2 , ' SUR TEN3=',AG3 

DO 200 J3=l , 5 

DO 190 J2»l,6 

DO 180 Jl=l,7 

AH-0. 5D0+(J3-1) *0 . 25D0 

AGRAV-GO* ( 1 . 0D0- (J2-1) *0 . 2 DO) 

AWN-0 . 25D0+ (Jl-1) *0. 5 DO 
C 

WRITE ( 14 , *) ' H* ' , AH , ' GRAV- ' , AGRAV 
WRITE (14,*) 'WAVE NUMBER-', AWN 
CALL DISP(ACFS) 

DO 50 J-1,5 
COEFF ( J) -ACFS ( J) 

50 CONTINUE 

C 

CALL DZPORC(NDEG, COEFF, ROOT) 

C 

WRITE ( 14 , *) 'RT1-' ,ROOT(l) 

WRITE ( 14 , *) ' RT2— # , ROOT ( 2 ) 

WRITE ( 14, *) 'RT3*' , ROOT ( 3 ) 

WRITE (14 , *) 'RT4- # , ROOT ( 4 ) 

180 CONTINUE 

190 CONTINUE 

200 CONTINUE 


CASE 3: AIR/SILICON OIL/AIR 
R2— DEN ( 3 ) 

R1«DEN(2) 

R3-DEN(3) 

AG2-GAM(3) 

WRITE (14,*) * • 

WRITE (14, *) # • 

WRITE ( 14, *) • AIR OVER SILICON OIL OVER AIR' 
WRITE (14,*) 'DEN2-' ,R2 
WRITE (14,*) 'DEN1-' ,R1 
WRITE (14,*) 'DEN3-' ,R3 

WRITE ( 14, *) 'SUR TEN2-' ,AG2, 'SUR TEN3-',AG3 

DO 300 K3-l,5 

DO 290 K2-l,6 

DO 280 Kl-1,7 

AH— 0 . 5D0+ (K3-1) *0.25D0 

AGRAV-GO* ( 1 . 0D0- (K2-1) *0 . 2D0) 
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c 


AWN=0 . 25D0+ (Kl-1) *0.5D0 


60 

C 

C 

280 

290 

300 


CASE 


C 


C 


C 


70 

C 


WRITE ( 14, *) ' H- ' , AH , ' GRAV= ' , AGRAV 
WRITE ( 14, *) 'WAVE NUMBER=',AWN 
CALL DISP(ACFS) 

DO 60 K-1,5 
CONTINUE 

CALL DZPORC (NDEG, COEFF, ROOT) 

WRITE ( 14 , * ) ' RT1— ' , ROOT ( 1 ) 

WRITE ( 14 , *) ' RT2— ' , ROOT ( 2 ) 

WRITE ( 14 , *) ' RT3- ' , ROOT ( 3 ) 

WRITE ( 14 , * ) 'RT4-' ,ROOT(4) 

CONTINUE 

CONTINUE 

CONTINUE 


: WATER/SILICON OIL/WATER 
R2-DEN(1) 

R1*DEN ( 2 ) 

R3-DEN(1) 

AG2-GAM(2) 

AG3«GAM(2) 

WRITE (14 , *) # ' 

WRITE (14,*)' ' 

WRITE(14,*) 'WATER OVER SILICON OIL OVER WATER' 
WRITE (14, *) 'DEN2-' ,R2 
WRITE (14,*) 'DEN1-',R1 
WRITE (14,*) 'DEN3-' ,R3 

DO 400 L3-l,5 
DO 390 L2-l,6 
DO 380 Ll-1,7 
AH-0.5D0+(L3-1)*0.25D0 
AGRAV-G0*(1.0D0-(L2-1) *0.2D0) 

AWN*0 . 2 5D0+ ( LI- 1 ) *0. 5D0 

WRITE (14, *) 'H-' , AH, ' GRAV- ', AGRAV 
WRITE (14,*) 'WAVE NUMBER-', AWN 
CALL DISP(ACFS) 

DO 70 L-1,5 
COEFF (L)-ACFS(L) 

CONTINUE 


177 



ooo ooo oo 


CALL DZPORC (NDEG, COEFF, ROOT) 

WRITE ( 14 , * ) ' RT1= ' , ROOT ( 1 ) 
WRITE (14, * ) 'RT2*' , ROOT ( 2 ) 
WRITE (14,*) 'RT3-' , ROOT (3) 
WRITE ( 14 , * ) 'RT4»' ,ROOT(4) 

380 CONTINUE 
390 CONTINUE 
400 CONTINUE 


CLOSE ( 14 ) 

STOP 

END 


SUBROUTINE DISP(ACFS) 

IMPLICIT REAL* 8 (A-H,0-Z) 

DIMENSION ACFS (5) 

COMMON/ DAT/ R 1 , R2 , R3 , AH, AG2 , AG3 , AGRAV, AWN 


COEFFICIENTS 

A-R1+R3 

B»( (AGRAV/AWN) *(R1-R3) ) - (AG3*AWN) 

C-R1+R2 

D-( (AGRAV/AWN) * (R2-R1) )-(AG2*AWN) 

W xo R2-Rl 

X»( (AGRAV/AWN) *(R2-R1) )-(AG2*AWN) 

Y»R1-R3 

Z-( (AGRAV/AWN) *(R3-R1) )+(AG3*AWN) 

ACFS ( 1) - (X+Z ) + (B*D*EXP ( 2 . 0D0*AWN*AH) ) 
ACFS(2)-0.0D0 

ACFS (3) *(W*Z) + (X*Y) + ( (A*D+B*C) *EXP(2 . 0D0*AWN*AH) ) 
ACFS ( 4 ) “0 . 0D0 

ACFS(5)-(W*Y)+(A*C*EXP(2. 0DO*AWN*AH) ) 

RETURN 

END 
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APPENDIX 5 


Generalized Eigenvalue Solution - Chap 3 

This program solves the large, sparse, generalized 
eigenvalue problem which is represented by equation (3.41) 
of Chapter 3 (periodic forcing case) . Truncation was made at 
N=|25| giving rise to an eigensystem of 204 equations. 

The complex eigenvalues are determined using ' DGVLCG ' of 
the IMSL library. The eigenvalues are the Floquet exponents of 
equations (3.35-3.38). Following computation of the 
eigenvalues, only the largest real component will be extracted 
for the data set. This represents the fastest growing Floquet 
exponent . 
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C PROGRAM: NEWPER FORTRAN 

C 

C 

C 3/21/90 
C 

C RE: PERIODIC MULTI-SLAB ANALYSIS 

C 

C 

C THIS PROGRAM IS STRUCTURED TO SOLVE A LARGE, SPARSE GENERALIZED 
C EIGEN VALUE MATRIX, RESULTING FROM AN ANALYSIS OF MULTI-LAYERED 
C SLABS OF LIQUID UNDER A NORMAL PERIODIC FORCING FUNCTION IN A 
C MICROGRAVITY ENVIRONMENT. 

C 

C FLOQUET THEORY WAS APPLIED, GENERATING A SYSTEM OF AN INFINITE 
C NUMBER OF LINEAR EQUATIONS. THIS SYSTEM WAS TRUNCATED AT 25 
C WHICH RESULTS IN 204 EQUATIONS. 

C 

C SINCE THE PROBLEM IS ESSENTIALLY A GENERALIZED COMPLEX 
C EIGENSYSTEM OF THE FORM, A*Z=W*B*Z, THE ROUTINE 'DGVLCG' OF THE 
C IMSL LIBRARY WILL BE CALLED TO SOLVE FOR THE EIGENVALUES. 

C THE EIGENVALUES (W) ARE THE FLOQUET EXPONENTS OF THE SYSTEM. 

C 

C FOLLOWING COMPUTATION OF THE EIGENVALUES, THE LARGEST POSITIVE 
C COMPONENT IS EXTRACTED FROM EACH ITERATION OF VARIOUS PARAMETER 
C THIS COMPONENT WILL DETERMINE THE STABILITY OF THE CONFIGURATIO 
C 

c 

c 

PARAMETER (N-204, NG-51) 

IMPLICIT REAL*8(A-H,0-Z) 

COMPLEX* 16 A(N,N) ,BLK(4,4) ,B(N,N) 

COMPLEX*16 ALPHA (N) , BETA (N) , EVAL(N) 

COMPLEX* 16 Y11,Y12, Y13, Y14 , Y21 , Y22 , Y23 , Y24, Y31, Y32, Y33, Y34 
COMPLEX* 16 Y41, Y42, Y43 , Y44, Yl, Y2, YB, YT, YDIA 
REAL RECO ( 4 ) ,Q(4) 

REAL AK 

EXTERNAL DGVLCG 
COMMON /WORKS P/ RWKSP 
REAL RWKSP( 332948) 

CALL IWKIN(332948) 

C 

LDA-N 

LDB-N 

C 

OPEN (UNIT- 16, STATUS “'NEW' , FILE-' FOR016 ' ) 

C 

C FILLING THE A AND B MATRICES WITH ZEROES PRIOR TO LOADING 
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C THE NON-ZERO TERMS 


DO 5 11*1 , N 
DO 4 JJ-1,N 
A(II,JJ)— (0. 0D0 , 0 . 0D0 ) 
B(II,JJ)=(O.ODO,O.ODO) 

4 CONTINUE 

5 CONTINUE 


THIS IS THE PARAMETER BLOCK. LOOPS ARE PERFORMED ON THE BOND 
NUMBERS ( B02 , B03 ) , THE FROUDE NUMBER (FR), AND THE WAVE NUMBER 
(K). VALUES OF THE DENSITY RATIOS (RH21,RH31) ARE SPECIFIED. 

RH21-10.0D0 
RH31— 0 . 001225D0 

RHDl-( DABS (RH2 1-1. 0D0 ) +DABS (RH3 1-1. 0D0) )/2.0D0 

DO 500 NJ«0,0 
DO 400 NT-0,0 
DO 300 IY-1,1 
WN-0 . 0D0 
DO 200 IX— 1 , 50 
WN-WN+0 . IDO 


FR-0. 01D0 

DO 160 NP-0,3 
NPP-NP+1 

B02— 10 . 0D0** ( -NP) 

B03-B02 

CALCULATING NON-ZERO ELEMENTS THAT WILL BE INSERTED INTO 
MATRICES A AND B. 

DO 20 M— 1 , NG 

SUB— ( (NG*1.0D0)+1.0D0)/2 .0D0 

CF- (M*l . 0D0) -SUB 
XI 3—1 . 0D0*WN*DEXP (WN) 

X14-1 . 0D0*WN*DEXP ( -1 . 0D0*WN) 

X23— 1 . 0D0*WN 
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X24=1.0D0*WN 

X31N=FR*WN*WN*DEXP(-1 . ODO*WN) * (RHD1 ) 

X3 1D=B02* (RH2 1+1 . ODO ) 

X31-X31N/X31D 

X34 - ( ( 1 . 0D0-RH2 1 ) *DEXP ( -2 . ODO*WN) ) / ( 1 . 0D0+RH2 1 ) 
X42=(-FR*WN*WN* (RHD1) ) / ( B03 * ( 1 . ODO+RH3 1 ) ) 

X43— ( 1 . ODO-RH3 1 ) / ( 1 . ODO+RH3 1) 

XB=- 1 . ODO* ( 1 . ODO-RH3 1 ) / ( 1 . ODO+RH31) 
XT=(DEXP(-2.0D0*WN) * (RH2 1-1 . ODO) ) / (RH2 1+1 . ODO ) 

XI- (FR* (1. ODO-RH21) )/ ( (1 . 0D0+RH21) *2 . ODO*DEXP(WN) ) 
X2-FR* ( 1 . 0D0-RH3 1)/(2.0D0*(1. ODO+RH3 1 ) ) 

X34M-X34*CF 

X43M-X43*CF 

ELEMENTS OF 4X4 SUBMATRIX FOR A GIVEN NG. 

XZER-0.0D0 

Y11«DCMPLX(XZER,CF) 

BLK(1,1)-Y11 

Y12*DCMPLX (XZER, XZER) 

BLK( I, 2 ) «Y12 

Y13-DCMPLX (X13 , XZER) 

BLK(1,3)-Y13 

Y14»DCMPLX(X14,XZER) 

BLK(1,4)-Y14 

Y21-DCMPLX (XZER, XZER) 

BLK(2 , 1) «Y21 

Y22*DCMPLX (XZER, CF) 

BLK(2,2)-Y22 

Y 2 3 “DCMPLX ( X2 3 , XZER) 

BLK(2,3)-Y23 

Y24-DCMPLX(X24 , XZER) 

BLK(2 , 4) "Y24 

Y3 1-DCMPLX (X3 1 , XZER) 

BLK(3,1)»Y31 

Y32-DCMPLX (XZER, XZER) 

BLK(3,2)-Y32 

Y3 3-DCMPLX ( XZER , CF) 

BLK(3,3)-Y33 

Y34-DCMPLX ( XZER , X3 4 ) 

BLK(3,4)«Y34 

Y4 1-DCMPLX ( XZER , XZER) 

BLK(4 , 1) -Y41 

Y42-DCMPLX (X42 , XZER) 

BLK(4 , 2) -Y42 

Y43-DCMPLX (XZER, X43 ) 

BLK(4 , 3) — Y43 

Y44-DCMPLX (XZER, CF) 
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BLK (4,4) =Y4 4 


LOADING NON-ZERO TERMS IN MATRIX A 

Yl-DCMPLX (XI , XZER) 

Y2=DCMPLX (X2 , XZER) 
NL»(4*M)-3 
NU«(4*M) 

Kl-0 

DO 15 I«NL,NU 

K2-0 

Kl-Kl+1 

DO 10 J«NL,NU 

K2-K2+1 

A(I, J) »BLK(K1, K2) 

IF(M.EQ. l.OR.M.EQ.NG) GO TO 8 
I F ( K1 . EQ . 3 ) THEN 
JB-I-6 
JF-I+2 
IP1-I+1 
JFP-JF+1 
JBP-JB+1 
A(I,JF)-Y1 
A(I,JB)-Y1 
A(IP1,JFP)-Y2 
A(IP1,JBP)-Y2 
ELSE 
END IF 
GO TO 10 

8 I F (M . EQ . 1 . AND . K1 . EQ . 3 ) THEN 

JF-I+2 
IP1-I+1 
JFP-JF+1 
A(I, JF)— Y1 
A(IP1, JFP) -Y2 
ELSE 
END IF 

IF(M.EQ.NG.AND.K1.EQ. 3) THEN 
JB-I-6 
IP1-I+1 
JBP-JB+1 
A(I,JB)-Y1 
A(IP1,JBP)-Y2 
ELSE 
END IF 
10 CONTINUE 
15 CONTINUE 
20 CONTINUE 
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LOADING NON-ZERO TERMS OF MATRIX B 


YB»DCMPLX(XB,XZER) 
YT-DCMPLX (XT, XZER) 
XN1— 1.0D0 

YDIA-DCMPLX (XN1 , XZER) 
NCT-0 

DO 30 L-1,NG 
DO 25 MOP-1,4 
NCT— NCT+1 
B (NCT, NCT) -YDIA 
IF(MOP.EQ. 4) THEN 
NM1-NCT-1 
B (NCT, NM1) -YB 
B(NM1,NCT) -YT 
ELSE 
END IF 
2 5 CONTINUE 
30 CONTINUE 


CALL DGVLCG (N, A, LDA, B, LDB, ALPHA, BETA) 

PROGRAM DGVLCG CALCULATES THE EIGENVALUES OF A GENERALIZED 
COMPLEX EIGENSYSTEM. 

THE EIGENVALUE (EVAL(N) ) IS COMPUTED BY DIVIDING COMPLEX 
VECTORS ALPHA (N) BY BETA(N) . 

THE EIGENVALUES ARE SWEPT OUT IN ORDER OF INCREASING SIZE OF 
THE REAL COMPONENT. THUS TO EXTRACT THAT VALUE, ONE NEEDS 
ONLY THE N-TH REAL VALUE OF EVAL. 

THIS LARGEST REAL COMPONENT (RECO) , IS THEN SENT TO A DATA 
FILE FOR VARIOUS PARAMETER VARIATIONS. 

DO 50 IM— 1 , N 

EVAL(IM) -ALPHA ( IM) /BETA ( IM) 

50 CONTINUE 

RECO (NPP) -EVAL (N) 


WRITE (16,2) WN , RECO (NPP) , B02 
2 FORMAT ( IX, F5 . 2 , E10 . 3 , F6 . 3 ) 

160 CONTINUE 
200 CONTINUE 
300 CONTINUE 
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400 CONTINUE 
500 CONTINUE 


CLOSE ( 16) 

STOP 

END 


185 



APPENDIX 6 


Standard Eigenvalue Problem - chap 3 

This program converts the generalized eigenvalue problem 
(of form $ X - A § X) to the standard form of g X » A X. This 
requires premultiplication of both sides by b” 1 using IMSL 
routine ' DLINGC' . 

The eigenvalues are calculated using routine * DEVLCG ' 
which utilizes a different algorithm than the generalized 
eigenvalue problem. The Floquet exponents as determined by 
both methods will be compared to check accuracy. 
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C PROGRAM 'LONG' 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


c 


THIS PROGRAM IS STRUCTURED TO SOLVE A LARGE. SPARSE GENERALIZED 
EIGEN VALUE MATRIX, RESULTING FROM AN ANALYSIS OF MULTI-LAYERED 
SLABS OF LIQUID UNDER A NORMAL PERIODIC FORCING FUNCTION IN A 
MICROGRAVITY ENVIRONMENT. 


THE PROBLEM CAN BE CONVERTED FROM A GENERALIZED TO A REGULAR 
EIGENVALUE PROBLEM BY PREMULTIPLYING BOTH SIDES BY BINV. THIS 
CAN BE CARRIED OUT BY IMPLEMENTING IMSL ROUTINES 'DLINCG' 
DMCRCR' , AND FINALLY THE EIGENVALUES CAN BE DETERMINED BY 
USING ROUTINE ' DEVLCG ' . 


PARAMETER (N-204, NG-51) 

IMPLICIT REAL* 8 (A-H, O-Z ) 

COMPLEX* 16 A(N, N) , BLK (4 , 4 ) , B (N, N) 

ALPHA (N) , BETA (N) , EVAL(N) 

Y11,Y12,Y13, Y14, Y21, Y22,Y23, Y24, Y31, Y3 2, Y3 3, Y3 4 
Y41,Y42,Y43,Y44,Y1,Y2,YB,YT, YDIA 
BINV (N, N) , C (N, N) 

EXTERNAL DGVLCG, DLINCG, DMCRCR, DEVLCG 
COMMON /WORKS P/ RWKSP 
REAL RWKSP (332948) 

CALL IWKIN (332948) 


COMPLEX*16 
COMPLEX* 16 
COMPLEX* 16 
COMPLEX* 16 


LDA-N 

LDB-N 

LDBINV-N 

NR-N 

LDC-N 


OPEN (UNIT-16, STATUS- 'NEW' , FILE-'FOR016' ) 
C 

C FILLING THE MATRICES WITH ZEROES 
C 

DO 5 II— 1,N 
DO 4 JJ— 1,N 
A (II, JJ) — (0 . 0D0, 0 . 0D0) 

B(II,JJ)— (0. 0D0 , 0 . 0D0) 

4 CONTINUE 

5 CONTINUE 
C 

C 

c 
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PARAMETER BLOCK 

BET=0 . 01D0 
WN-1. 0D0 
DO 300 IY«1,1 
WN-WN+0. 5 DO 
DO 200 IX-1,1 
BET-BET* (5. ODO*IX) 

T2=0 . 001D0 
T3-0.001D0 
RH2 1-0 . 00 12 2 5 DO 
RH31— 0 . 001225D0 

RHD1- (DABS (RH21-1. 0D0) +DABS (RH31-1. 0D0) )/2 . 0D0 


CALCULATING NON-ZERO ELEMENTS 
DO 20 M»1,NG 

SUB- ( (NG*1 . 0D0) +1 . 0D0) /2 . 0D0 

CF— (M*l . 0D0) -SUB 

X13— 1 . 0D0*WN*BET*DEXP (WN) 

X14-1 . 0D0*WN*BET*DEXP (-1 . 0D0*WN) 

X23--1 . 0D0*WN*BET 
X24-1 . 0D0*WN*BET 

X31N— WN*WN*DEXP(-1. 0D0*WN) * (RHD1) 

X31D-T2* (RH21+1 . 0D0) 

X31-X31N/X31D 

X34— ( -1 . 0D0*DEXP(-1 . 0D0*WN) * (RH21-1 . 0D0) *CF) / (RH2 1+1 . 0D0 ) 
X42- ( -1 . 0D0*WN*WN* (RHD1 ) ) / (T3 * ( 1 . 0D0+RH3 1 ) ) 

X43— ( ( 1 . 0D0-RH31) *CF)/ ( 1 . 0D0+RH31) 

XB— 1 . 0D0* ( 1 . 0D0-RH3 1) / ( 1 . 0D0+RH3 1) 

XT-(DEXP(-1.0D0*WN) * (RH2 1-1 . 0D0) ) / (RH21+1 . 0D0) 

XI— 1.0D0*(RH21-1. 0D0)/ (2. 0D0* (RH21+1. 0D0) *DEXP(WN) ) 

X2— (1. 0D0-RH31)/ (2 . 0D0* ( 1 . 0D0+RH31) ) 

X34M-X34*CF 

X43M-X43*CF 

ELEMENTS OF 4X4 SUBMATRIX 

XZER-0 . 0D0 

Yll— DCMPLX (XZER, CF) 

BLK( 1 , 1) — Yll 

Y12-DCMPLX (XZER, XZER) 

BLK( 1,2) — Y12 

Y13-DCMPLX (X13 , XZER) 

BLK( 1 , 3 ) -Y13 

Y14-DCMPLX (X14 , XZER) 


188 



ooo 


BLK(1,4)«Y14 

Y21-DCMPLX (XZER, XZER) 
BLK(2,1)«Y21 

Y22=DCMPLX (XZER, CF) 
BLK(2,2)-Y22 

Y23*DCMPLX(X23 , XZER) 
BLK(2,3)-Y23 

Y24-DCMPLX (X24 , XZER) 
BLK(2,4)-Y24 

Y3 1-DCMPLX (X3 1 , XZER) 
BLK(3,1)-Y31 

Y32-DCMPLX (XZER, XZER) 

BLK(3 , 2) *Y32 

Y33»DCMPLX(XZER,CF) 

BLK(3,3)-Y33 

Y34“DCMPLX (XZER, X34 ) 
BLK(3,4)-Y34 

Y41-DCMPLX( XZER, XZER) 
BLK(4,1)-Y41 

Y42-DCMPLX (X42 , XZER) 
BLK(4,2)-Y42 

Y43*DCMPLX (XZER, X4 3 ) 
BLK(4,3)-Y43 

Y44-DCMPLX(XZER,CF) 

BLK(4,4)-Y44 

LOADING NON-ZERO TERMS IN MATRIX A 

Y1“DCMPLX (XI , XZER) 

Y2“DCMPLX ( X2 , XZER) 
NL-(4*M)-3 
NU-(4*M) 

Kl-0 

DO 15 I«NL,NU 

K2-0 

Kl-Kl+1 

DO 10 J*NL,NU 

K2-K2+1 

A(I,J)-BLK(K1,K2) 
IF(M.EQ.l.OR.M.EQ.NG) GO TO 8 
IF(K1.EQ.3)THEN 
JB-I-6 
JF-I+2 
IP1-I+1 
JFP-JF+1 
JBP-JB+1 
A(I,JF)«Y1 
A(I,JB)-Y1 
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A ( IP1 , JFP) =Y2 
A ( IP1 , JBP) =*Y2 
ELSE 
END IF 
GO TO 10 

IF(M. EQ. 1 . AND. K1 . EQ. 3 ) THEN 
JF-I+2 
IP1-I+1 
JFP-JF+l 
A ( I , JF) =Y 1 
A ( IP1 , JFP) =Y2 
ELSE 
END IF 

IF (M. EQ.NG. AND. K1 . EQ. 3 ) THEN 
JB-I-6 
IP1-I+1 
JBP-JB+1 
A(I,JB)-Y1 
A ( IP1 , JBP) — Y2 
ELSE 
END IF 
10 CONTINUE 
15 CONTINUE 
20 CONTINUE 


LOADING NON-ZERO TERMS OF MATRIX B 

Y B- DCMPLX ( XB , XZ ER ) 

YT-DCMPLX (XT, XZER) 

XN1— 1.0D0 

YDIA-DCMPLX ( XN1 , XZER) 

NCT-0 

DO 30 L-1,NG 
DO 25 MOP-1,4 
NCT-NCT+1 
B(NCT,NCT) -YDIA 
IF (MOP. EQ. 4) THEN 
NM1-NCT-1 
B(NCT,NM1)-YB 
B(NM1,NCT)-YT 
ELSE 
END IF 
25 CONTINUE 
30 CONTINUE 

COMPUTING THE INVERSE OF B 
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c 


CALL DLINCG (N, B, LDB , BINV, LDBINV) 


MULTIPLYING BINV AND A 

CALL DMCRCR (NR, NR, BINV, LDA, NR, NR, A, LDB, NR, NR, C, LDC) 

SOLVING FOR THE EIGENVALUES 

CALL DEVLCG (N, C, LDC, EVAL) 

DO 50 IM»1,N 

WRITE ( 16, *) ' EVAL- ' , EVAL ( IM) 

50 CONTINUE 

00 CONTINUE 
00 CONTINUE 


CLOSE (16) 

STOP 

END 
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APPENDIX 7 


Determinant Calculation 


This program takes the eigenvalues which were determined 
by Appendix 5 and substitutes them into equation (3.39). The 
determinant of the resulting matrix is calculated to check for 
accuracy of the eigenvalues. The determinant should equal zero 
if the eigenvalues are accurate. 

Routine 'DLFTCG' is utilized for LU factorization of the 
matrix and computation of the determinant. 
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C THIS PROGRAM ' DET FORTRAN' WILL CALCULATE THE DETERMINANT OF 

THE SPARSE MATRIX TO DETERMINE THE ACCURACY OF THE SOLUTION. 

THE ROUTINE ' DLFDCG ' OF THE IMSL PACKAGE WILL BE UTILIZED. THE 
LU FACTORIZATION WILL BE NEEDED AND WILL PROVIDED BY CALLING THE 
' DLFTCG ' ROUTINE. 


THIS PROGRAM IS STRUCTURED TO SOLVE A LARGE, SPARSE GENERALIZED 
EIGEN VALUE MATRIX, RESULTING FROM AN ANALYSIS OF MULTI-LAYERED 
SLABS OF LIQUID UNDER A NORMAL PERIODIC FORCING FUNCTION IN A 
MICROGRAVITY ENVIRONMENT. 

SINCE THE PROBLEM IS ESSENTIALLY A GENERALIZED COMPLEX 
EIGENSYSTEM OF THE FORM, A*Z=W*B*Z, THE ROUTINE ' DGVLCG' OF THE 
IMSL LIBRARY WILL BE CALLED TO SOLVE FOR THE EIGENVALUES. 


PARAMETER (N-204, NG-51) 

IMPLICIT REAL*8 ( A-H , O-Z ) 

COMPLEX* 16 A(N,N) ,BLK(4, 4) ,B(N,N) 

COMPLEX* 16 ALPHA (N) , BETA (N) , EVAL(N) 

COMPLEX*16 FAC(N,N) , DET1 , AMAT (N, N) 

COMPLEX* 16 Y11,Y12,Y13,Y14,Y21,Y22,Y23,Y24 , Y31, Y32, Y33 , Y34 
COMPLEX* 16 Y41, Y42 , Y43 , Y44 , Yl, Y2 , YB, YT, YDIA 
REAL* 8 DET2 
INTEGER IPVT(N) 

EXTERNAL DGVLCG, DLFDCG, DLFTCG 
COMMON /WORKS P/ RWKSP 
REAL RWKSP (332948) 

CALL IWKIN (332948) 

LDA-N 

LDB-N 

LDFAON 

OPEN (UNIT-19 , STATUS- 'NEW ' , FILE- ' FOR019 ' ) 

FILLING THE MATRICES WITH ZEROES 

DO 5 II— 1,N 
DO 4 JJ-1 , N 
A (II, JJ)-(0.0D0,0.0D0) 

B( II , JJ) — (0 . 0D0 , 0 . 0D0) 

4 CONTINUE 

5 CONTINUE 
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PARAMETER BLOCK 

BET-0. 01D0 

WN-1.0D0 

DO 300 IY-1,1 

WN-WN+0.5D0 

DO 200 IX-1,1 

BET-BET* (5. 0D0*IX) 

T2-0.001D0 
T3-0.001D0 
RH21-0 . 001225DO 
RH31-0 . 001225D0 

RHD1— (DABS (RH21-1 . 0D0) + DABS (RH31-1.0D0) )/2 . 0D0 


CALCULATING NON-ZERO ELEMENTS 
DO 20 M— 1,NG 

SUB— ( (NG*1 . 0D0) +1 . 0D0) /2 . 0D0 

CF— (M*l . 0D0) -SUB 

X13--1 . 0D0*WN*BET*DEXP (WN) 

X14-1 . 0D0*WN*BET*DEXP ( -1 . 0D0*WN) 

X23— 1.0D0*WN*BET 
X24— 1 . 0D0*WN*BET 

X3 1N=WN*WN*DEXP(-1 . 0D0*WN) * (RHD1) 

X31D-T2* (RH21+1 . 0D0) 

X31-X31N/X31D 

X34— (-1 . 0D0*DEXP(-1 . 0D0*WN) * (RH21-1 . 0D0) *CF)/ (RH21+1 . 0D0) 
X42- (-1 . 0D0*WN*WN* (RHD1) )/ (T3* (1. 0D0+RH31) ) 

X43- ( ( 1 . 0D0-RH31) *CF) / ( 1 . 0D0+RH31) 

XB— 1 . 0D0* ( 1 . 0D0-RH3 1) / ( 1 . 0D0+RH3 1 ) 

XT— (DEXP(-1. 0D0*WN) * (RH21-1 . 0D0) )/ (RH21+1 . 0D0) 

XI— 1.0D0*(RH21-1.0D0)/ ( 2 . 0D0* (RH21+1 . 0D0) *DEXP(WN) ) 

X2— ( 1 . 0D0-RH3 1) / ( 2 . 0D0* ( 1 . 0D0+RH3 1) ) 

X34M— X34*CF 
X43M-X43*CF 

ELEMENTS OF 4X4 SUBMATRIX 

XZER-0.0D0 
Y 1 1-DCMPLX ( XZER , CF ) 

BLK(1, 1) — YI1 

Y12-DCMPLX( XZER, XZER) 

BLK( 1 , 2 ) -Y12 
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Y13-DCMPLX(X13 ,XZER) 

BLK(1, 3)-Y13 

Y14— DCMPLX (X14 , XZER) 
BLK(1,4)»Y14 

Y21=DCMPLX( XZER, XZER) 

BLK(2 , 1) -Y21 

Y2 2— DCMPLX (XZER, CF) 
BLK(2,2)«Y22 

Y23«DCMPLX(X23 , XZER) 

BLK(2, 3)-Y23 

Y 2 4 -DCMPLX ( X 2 4 , X Z ER ) 
BLK(2,4)-Y24 

Y3 1-DCMPLX (X3 1 , XZER) 
BLK(3,1)-Y31 

Y32-DCMPLX (XZER, XZER) 
BLK(3,2)-Y32 

Y33-DCMPLX (XZER, CF) 
BLK(3,3)-Y33 

Y34-DCMPLX(XZER,X34) 

BLK( 3 , 4 ) -Y34 

Y41-DCMPLX (XZER, XZER) 

BLK(4 , 1)-Y41 

Y42-DCMPLX (X42 , XZER) 
BLK(4,2)-Y42 

Y4 3 -DCMPLX (XZER, X43 ) 
BLK(4,3)«Y43 

Y4 4 -DCMPLX (XZER, CF) 

BLK (4,4) — Y44 

LOADING NON-ZERO TERMS IN MATRIX A 

Y1-DCMPLX(X1,XZER) 

Y2-DCMPLX (X2 , XZER) 
NL»(4*M)-3 
NU- (4*M) 

Kl-0 

DO 15 I-NL, NU 

K2-0 

Kl-Kl+1 

DO 10 J-NL, NU 

K2-K2+1 

A(I, J) — BLK(K1 , K2) 

IF(M.EQ. l.OR.M.EQ.NG) GO TO 8 
IF(K1.EQ.3)THEN 
JB-I-6 
JF-I+2 
IP1-I+1 
JFP-JF+1 
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JBP-JB+l 


A(I, JF)-Y1 
A(I, JB)»Y1 
A ( IP1 , JFP) =Y2 
A(IP1, JBP)*Y2 
ELSE 
END IF 
GO TO 10 

8 IF(M.EQ. 1.AND.K1.EQ. 3) THEN 
JF-I+2 
IP1-I+1 
JFP-JF+1 
A(I, JF)«Y1 
A(IP1 # JFPJ-Y2 
ELSE 
END IF 

IF (M . EQ . NG . AND . K1 . EQ . 3 ) THEN 
JB-I-6 
IP1-I+1 
JBP-JB+l 
A(I , JB) — Y1 
A(IP1 , JBP) — Y2 
ELSE 
END IF 
10 CONTINUE 
15 CONTINUE 
20 CONTINUE 


LOADING NON-ZERO TERMS OF MATRIX B 

YB-DCMPLX (XB, XZER) 

YT-DCMPLX (XT, XZER) 

XN1— 1.0D0 

YDIA-DCMPLX (XN1 , XZER) 

NCT— 0 

DO 30 L— 1 ,NG 
DO 25 MOP-1,4 
NCT— NCT+1 
B(NCT,NCT)-YDIA 
IF (MOP. EQ. 4) THEN 
NM1-NCT-1 
B(NCT,NM1)-YB 
B(NM1,NCT) -YT 
ELSE 
END IF 

25 CONTINUE 
30 CONTINUE 
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CALL DGVLCG (N, A, LDA , B, LDB , ALPHA, BETA) 
DO 50 IM*1,N 

EVAL ( IM) -ALPHA ( IM) /BETA ( IM) 

50 CONTINUE 


200 CONTINUE 
300 CONTINUE 


COMPUTE NEW MATRIX A=A-W*B 

DO 500 MM»1,N 
DO 400 IL-1,N 
DO 350 JL-1,N 
AMAT ( IL, JL) -A ( IL, JL) 

A(IL, JL) -AMAT (IL,JL) -(EVAL (MM) *B(IL,JL) ) 
350 CONTINUE 
400 CONTINUE 

FACTORING MATRIX A 

CALL DLFTCG (N, A, LDA, FAC, LDFAC, I PVT) 


COMPUTE THE DETERMINANT OF THE FACTORED MATRIX 

CALL DLFDCG (N, FAC, LDFAC, IPVT, DET1, DET2 ) 


WRITE ( 19 , * ) DET1 , DET2 

DO 460 IR— 1 , M 
DO 450 JR-1 , N 
A ( IR , JR) -AMAT ( IR , JR) 
450 CONTINUE 
460 CONTINUE 
500 CONTINUE 
CLOSE (19) 

STOP 

END 
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One Interface Solution -Chap 3 

As discussed in section 3.4a, a limit approximation is 
compared to the 1 interface system of equations (see Appendix 
2). The two linear equations (A2.8,A2.9) are solved using 
Floquet analysis resulting in a standard eigensystera. 

Routine ' DEVLCG ' is used to determine the eigenvalues of 
the 1 interface configuration. The results are used to compare 
the limit approximation for the 2 interface system. 


198 



oo ooo oo o oo ooooooooooo 


C THIS IS A CHECK FOR ONE INTERFACE 


THIS PROGRAM SOLVES THE EIGENSYSTEM RESULTING FROM A ONE 
INTERFACE CONFIGURATION. 

THE SYSTEM REDUCES TO TWO LINEAR DIFFERENTIAL EQUATIONS WHICH 
ARE SOLVED VIA FLOQUET ANALYSIS IN THE SAME MANNER AS THE 
TWO INTERFACE SYSTEM. 

A STANDARD EIGENVALUE PROBLEM IS OBTAINED AND IS SOLVED USING 
' DEVLCG ' FROM THE IMSL LIBRARY. 


PARAMETER (N-102 , NG-5 1 ) 

IMPLICIT REAL* 8 (A-H,0-Z) 

COMPLEX* 16 A(N,N) , BLK(2,2) 

COMPLEX* 16 EVAL(N) 

COMPLEX* 16 Yll, Y12 , Y21, Y22, YOUT 
EXTERNAL DEVLCG 
COMMON /WORKS P/ RWKSP 
REAL RWKSP (3 32948) 

CALL IWKIN(332948) 


LDA-N 

OPEN (UNIT-16, STATUS- ' NEW ' , FILE- ' FORO 16 ' ) 

FILLING THE A MATRIX WITH ZEROES 
DO 5 1-1 , N 
DO 4 J-1,N 
A(I , J) — (0 . 0D0, 0 . 0D0) 

4 CONTINUE 

5 CONTINUE 


PARAMETER BLOCK 

RH21-0. 001225D0 
DO 500 NJ-0,0 
Q«1.0D0*10.0D0**NJ 
WN— 0 . 0D0 
DO 400 NT— 1, 50 
WN-WN+0 . 05D0 

CALCULATING NON-ZERO TERMS OF A 
DO 100 M-l , NG 

SUB— ( (NG*1 . 0D0) +1 . 0D0) /2 . 0D0 
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CF-(M*1. 0D0) -SUB 
X21— Q*WN*WN 
X12-WN 

XOUT» ( 1 . 0D0-RH2 1 ) / ( 2 . 0D0* ( 1 . 0D0+RH21) ) 
FC=-CF 


ELEMENTS OF 2X2 BLOCK 
XZER-O.ODO 
Yll-DCMPLX (XZER, FC) 
BLK( 1 , 1) =Y11 
Y12»DCMPLX (X12 , XZER) 
BLK(1,2)-Y12 
Y2 1-DCMPLX (X2 1 , XZER) 
BLK(2, 1) »Y21 
Y22-DCMPLX(XZER,FC) 
BLK(2 , 2 ) *Y22 
YOUT-DCMPLX (XOUT, XZER) 

LOADING TERMS OF A 
MP»(2*M)-1 
MP1-MP+1 
MPM-MP1-3 
MPP-MP1+1 
A (MP, MP) "BLK (1,1) 

A (MP, MP1) ■BLK( 1 , 2 ) 
A(MP1,MP)-BLK(2,1) 

A (MP1 , MP1 ) «BLK (2,2) 
IF(M.NE. 1)THEN 
A (MP1 , MPM) “YOUT 
ELSE 
END IF 

IF(M.NE.NG)THEN 
A (MP1 , MPP) "YOUT 
ELSE 
END IF 

100 CONTINUE 


CALL DEVLCG (N, A, LDA, EVAL) 


C 

C 


RECO-EVAL(N) 

UNUM-1.0 

WN-WN+0.05D0 
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WRITE ( 1 6 , 2 8 ) WN , UNUM , RECO 

28 FORMAT (1X,F5.2,F6.3,E10.3) 
C 

WN=WN-0.05D0 

C 

400 CONTINUE 
500 CONTINUE 
C 

CLOSE (16) 

STOP 

END 
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Time Response of Interfaces - Chap 4 

For the non-periodic forcing case, a system of linear 
differential equations in terms of (E,F,E,F) is obtained 
(equations (4.22,4.23)). E and F are the displacements of the 
upper and lover interfaces, respectively. 

The system of equations is integrated numerically using 
' DIVPAG ' of the IMSL library. This routine utilizes Gear's to 
solve for the time-dependent coefficients. 
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c 


THIS FILE IS CSLAB FORTRAN 


C THIS PROGRAM INTEGRATES THE FOURTH ORDER MULTI -SLAB 
C FORCED LINEAR SYSTEM, CONSISTING OF E ( T) , F , DE/ DT , DF/ DT 
C ( ALL FCNS . OF T) 

C THE FORCING IS **NOT*** PERIODIC 

C GEAR'S METHOD IS USED — PROB MAY BE STIFF. 

C IMSL LIBRARIES ARE USED. 

C PARAMETER VALUES ARE CHOSEN. 

C 

C 

INTEGER NEQ , NPARAM 
PARAMETER (NPARAM-50 , NEQ-4 ) 

C 

INTEGER IDO, IEND, IMETH, INORM,NOUT 

REAL*8 A ( 1 , 1) ,FCN,FCNJ,HINIT,PARAM (NPARAM) , TOL,T, TEND, Y (NEQ) 
REAL*8 AK, B2,B3,R21,R31, RHD1 
REAL* 8 EE ( 3 ) ,FF(3) 

EXTERNAL FCN, DIVPAG , SSET, UMACH 
C 

COMMON/DAT/ AK , B2 , B3 , R2 1 , R3 1 , RHD1 


OPEN (UNIT-2 4, STATUS- 'NEW' , FILE-'FOR024 ' ) 

COUNTER 

KK-0 

JEDE-0 

PARAMETER VALUES 
AK-1.0D0 
B2-1.0D0 
B3-1.0D0 
R21-0 . 001225D0 
R31-0. 001225D0 

RHD1- (DABS (R21-1 . 0D0) +DABS (R3 1-1 . 0D0) ) /2 . ODO 


DO 200 ID-0,1 

B2— 10. ODO** (-ID) 
B3-B2 

PP-ID+1 


SET 


C 


INITIAL CONDITIONS 
T-0.00D0 
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E=Y1, DE/DT-Y2, F-Y3 , DF/DT=»Y4 

CC— 0.05D0 
DD-0.05D0 

Y ( 1) *0 . 00D0 
Y ( 2 ) — AK*CC*DEXP ( -AK) 

Y ( 3 ) =0 . 00D0 
Y (4 ) =AK*DD 

SET PROGRAM SWITCHES/VALUES 
HINIT-0.0001D0 
I NORM” 1 
IMETH-2 

CALL SSET(NPARAM, 0.0, PARAM, 1) 
PARAM ( 1 ) — HINIT 
PARAM ( 10 )-INORM 
PARAM ( 12 )«IMETH 
SET ERROR TOLERANCE 
TOL-l.OD-5 


WN-AK 

B02-B2 

B03-B3 

IDO-1 

DO 100 11-1,5000 
TEND— 0 . 01D0*II 

CALL DIVPAG ( IDO , NEQ , FCN , FCNJ , A , T , TEND , TOL , PARAM , Y ) 
KK-KK+1 

IF (KK.EQ.10) GO TO 50 
GO TO 60 
50 CONTINUE 


JEDE-JEDE+1 

TT-T 

GRAV-G(T) 

IF ( JEDE.GE. 150) GRAV-0 . 0D0 
WRITE (24,28)T,Y(1) ,Y(3) , B02 
28 FORMAT ( IX , 4 F9 . 4 ) 

C 

KK-0 

60 CONTINUE 
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100 

c 

c 


c 

200 

C 


C 


c 


c 

c 


c 

c 


CONTINUE 


RELEASE WORKSPACE 
IDO=3 

CALL DIVPAG ( IDO, NEQ, FCN , FCNJ , A, T, TEND , TOL, PARAM, Y) 
CONTINUE 

CLOSE (24) 

END 


SUBROUTINE FCN (NEQ, T, Y , YP) 
INTEGER NEQ 


REAL* 8 T, Y(NEQ) ,YP(NEQ) 

REAL* 8 CF1N, CF1D, CF1 , CF2N , CF2D, CF2 , CBLK2 , CBLK3 , CFA, CFB 
REAL* 8 AK, R2 , R3 , RHD1 , B2 , B3 
COMMON/ DAT/ AK, B2 , B3 , R2 , R3 , RHD1 


CF1N— 2 . 0 DO * AK* DEXP ( -AK) 

CF1D- ( 1 . 0D0+R3 ) * ( 1 . 0D0+R2 ) + ( 1 . 0D0-R3 ) * (R2-1 . 0D0 ) 
& *DEXP ( -2 . 0D0*AK) 

CF1-CF1N/CF1D 

CF2N-AK*( (R2-1.0D0) *DEXP(-2 . 0D0*AK) - (R2+1 . ODO) ) 

CF2D-CF1D 

CF2-CF2N/CF2D 

CBLK2— (AK*AK/B2) * (RHD1) - (R2-1 . ODO) *G(T) 
CBLK3-(AK*AK/B3) * (RHD1 ) - ( 1. 0D0-R3) *G(T) 

CFA- ( 1 . 0D0-R3 ) / ( 1 . 0D0+R2 ) 

CFB-AK/ ( 1 . 0D0+R2 ) 


YP( 1) — Y (2) 

YP(2) — (CF1*CFA*DEXP(-AK) -CFB) *CBLK2*Y(1) 

& +CF1*CBLK3*Y (3) 

YP(3) — Y (4) 

YP(4)— (DEXP(-AK) * (CF2*CFA-CFB) *CBLK2) *Y ( 1) +CF2*CBLK3*Y ( 3 ) 

RETURN 

END 


FUNCTION G(T) 

REAL* 8 G,T 

IF(T.LT.l.O) G-0.0D0 

IF (T.GE. 1 . 0 . AND. T. LT. 2 . 0) G— T-l. ODO 

IF (T.GE. 2 . 0 . AND.T. LT. 3 . 0) G-l . ODO 

IF(T .GE. 3 . 0. AND.T. LT. 5 . 0) G--T+4 . ODO 

IF (T. GE . 5 . 0 . AND. T. LT. 6 . 0) G=-l . ODO 
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IF (T. GE . 6 . 0 . AND. T. LT . 7 . 0) G=T-7 . 0D0 

IF(T.GE.7 . 0D0) G=0 . 0D0 
RETURN 
END 
C 

SUBROUTINE FCNJ (NEQ, T, Y , YP) 

C DUMMY 

RETURN 
END 
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